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Abstract 

In  marine  ecology,  the  variability  of  the  physical  environment  is  often  considered  a  main 
determinant  of  biological  pattern.  A  common  approach  to  identifying  key  environmental 
forcings  is  to  match  scales  of  variability:  fluctuations  of  a  biological  variable  at  a  particular 
frequency  are  attributed  to  forcing  by  the  physical  environment  at  a  similar  frequency. 
In  nonlinear  systems,  however,  different  scales  of  variability  interact  and  forcing  at  one 
frequency  can  produce  variability  at  a  different  frequency. 

The  general  theme  of  this  dissertation  regards  the  interplay  of  scales  in  nonlinear  ecolog¬ 
ical  systems,  with  an  emphasis  on  the  mismatch  of  scales  between  biological  variables  and 
environmental  forcings  in  the  plankton.  The  approach  is  theoretical:  I  use  simple  models 
to  identify  conditions  leading  to  such  a  mismatch.  The  models  are  motivated  by  plank¬ 
tonic  systems  and  focus  on  one  ubiquitous  nonlinear  ecological  interaction,  that  between  a 
consumer  and  a  resource. 

This  work  is  organized  in  three  main  parts  as  follows.  In  the  first  part,  I  consider  the 
interaction  between  a  phytoplankton  population  and  a  limiting  nutrient  resource.  Most 
models  for  this  interaction  consider  all  cells  as  equal  and  group  them  under  a  single  vari¬ 
able,  the  total  biomass  or  cell  density.  They  do  not  take  into  account  any  population 
heterogeneity  resulting  from  the  life  histories  of  individual  cells.  However,  single  cells  do 
have  life  histories:  each  cell  progresses  through  a  determinate  sequence  of  events  preceding 
cell  division  and  the  population  is  distributed  in  stages  of  the  cell  cycle.  I  incorporate  this 
distribution  (i.e.  population  structure),  as  well  as  observations  on  resource  control  of  cell 
cycle  progression,  into  chemostat  models  for  the  phytoplankton-nutrient  interaction.  Simu¬ 
lation  results  demonstrate  that  the  population  structure  can  generate  oscillatory  dynamics 
under  a  constant  nutrient  supply,  and  that  such  oscillations  are  important  to  population 
dynamics  under  a  variable  nutrient  supply.  Specifically,  for  a  periodic  resource  supply,  the 
population  displays  an  aperiodic  response  with  frequencies  different  from  that  of  the  forcing. 
I  then  show  that  a  chemostat  model  without  population  structure  (the  Droop  equations) 
does  not  exhibit  this  transfer  of  variability:  a  periodic  nutrient  supply  produces  a  periodic 
population  response  of  exactly  the  same  frequency. 

In  the  second  part,  I  consider  a  predator  and  a  prey  that  interact  and  diffuse  along  an 
environmental  gradient.  The  model  is  a  reaction-diffusion  equation,  a  type  of  model  used  in 
biological  oceanography  for  planktonic  interactions  in  turbulent  flows.  I  demonstrate  that 
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weak  diffusion  along  a  spatial  gradient  may  drive  an  otherwise  periodic  system  into  complex 
temporal  dynamics  that  include  chaotic  behavior.  I  provide  evidence  for  a  quasiperiodic 
route  to  chaos  as  the  diffusion  rate  decreases.  Then,  I  focus  on  the  spatial  properties  of 
the  gradient  and  their  consequences  for  the  spatio-temporal  dynamics  of  the  system.  In 
particular,  I  ask:  how  do  the  spatial  patterns  of  the  populations  compare  to  the  underlying 
environmental  gradient  in  the  different  dynamic  regimes  (periodicity,  quasiperiodicity,  and 
chaos)?  I  show  that  the  spatial  patterns  of  predator  and  prey  can  differ  strongly  from  the 
environmental  gradient.  In  the  route  to  chaos,  as  diffusion  becomes  weaker,  this  difference 
is  magnified  and  the  populations  display  smaller  spatial  scales. 

In  the  work  summarized  so  far,  nonlinearity  leads  to  variability  in  biological  variables 
at  scales  not  present  in  the  environmental  forcings.  In  the  third  part  of  this  work,  I  con¬ 
sider  another  consequence  of  the  transfer  of  variability  in  nonlinear  systems:  the  lack  of  a 
dominant  scale.  Patterns  that  lack  a  dominant  scale  but  exhibit  scale  similarity  are  known 
as  fractals.  The  characterization  of  numerical  quantities  that  vary  intermittently  has  mo¬ 
tivated  a  generalization  of  fractals  known  as  multifractals.  Here,  I  give  a  first  application 
of  multifractals  to  biological  oceanography.  I  analyze  an  acoustic  data  set  on  zooplank¬ 
ton  biomass  to  describe  the  distribution  in  time  of  the  total  variability  in  the  data.  This 
distribution  is  highly  intermittent:  extreme  localized  contributions  account  for  a  large  pro¬ 
portion  of  total  variability.  I  show  that  multifractals  provide  a  good  characterization  of 
such  variability. 

Dissertation  Advisor:  Hal  Caswell 
Title:  Senior  Scientist,  WHOI 
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Chapter  1 


Introduction 


The  relations  between  the  two  of  them  must  have  been  fascinating. 

For  things  are  not  what  they  seem. 

— Saul  Bellow.  It  All  Adds  Up.  From  the  Dim  Past  to  the  Uncertain  Future. 

Nonlinearity,  which  is  common  in  ecological  interactions,  creates  a  rich  array  of  possible 
dynamics.  Of  these,  chaos  occupies  the  center  stage.  First  encountered  at  the  turn  of  this 
century  (Duhem,  1906;  Hadamard,  1898;  Poincare,  1908),  chaos  raised  the  unexpected  and 
somewhat  disturbing  possibility  of  aperiodic  dynamics  with  sensitivity  to  initial  conditions 
in  deterministic  systems.  Many  decades  passed  before  its  definite  rediscovery  in  the  work 
of  Lorenz  (1963),  Ruelle  and  Takens  (1971  ),  and  May  (1974).  A  watershed  followed  in 
nonlinear  dynamics  research.  Hassell  et  al.  (1976)  and  Schaffer  and  Kot  (1985,  1986) 
first  applied  approaches  from  nonlinear  dynamical  systems  to  ecological  data;  the  debate 
continues  on  the  relevance  of  chaos  to  natural  systems  (  Hastings  et  al.,  1993;  Ellner  and 
Turchin,  1995).  Chaos  remains  a  fascinating  concept:  it  introduces  the  notions  of  strange 
attractors  with  fractal  geometries,  complex  dynamics  in  simple  systems,  and  sensitivity  to 
initial  conditions  in  spite  of  determinism. 

It  is,  however,  on  a  different  but  related  property  of  nonlinearity  that  I  wish  to  focus 
my  attention  here.  This  property  involves  the  concept  of  scale,  which  is  used  throughout 
this  work  to  denote  the  dominant  period  (i.e.  the  inverse  of  the  dominant  frequency) 
or  the  dominant  wavelength  in  the  variance  of  a  temporal  or  spatial  quantity  of  interest. 
In  nonlinear  dynamical  systems,  different  scales  of  variability  interact  and  forcing  at  one 
scale  can  produce  variability  at  a  different  scale.  An  example  of  this  transfer  of  variability 
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between  scales  is  found  in  the  response  of  nonlinear  systems  to  forcing  by  a  single  temporal 
frequency.  When  this  response  is  chaotic,  it  exhibits  all  frequencies  (i.e,  a  continuous  power 
spectrum)  in  spite  of  the  the  single  forcing  frequency. 

How  does  the  interplay  of  scales  in  nonlinear  systems  matter  to  ecology,  and  more  specif¬ 
ically  to  plankton  dynamics?  Ecology  is  concerned  with  the  spatial  and  temporal  patterns 
of  populations  and  communities,  and  with  identifying  the  processes  responsible  for  such  pat¬ 
terns.  In  marine  ecology,  the  variability  of  the  physical  environment  is  often  considered  a 
main  determinant  of  biological  pattern  (Steele  and  Henderson,  1994).  A  common  approach 
to  identifying  key  environmental  forcings  is  to  match  scales  of  variability:  fluctuations  of  a 
biological  variable  at  a  particular  frequency  are  attributed  to  forcing  by  the  physical  environ¬ 
ment  at  a  similar  frequency.  This  approach  has  been  used  extensively  and  often  succesfully 
in  the  study  of  planktonic  systems  (Denman,  1994;  Denman  and  Powell,  1984),  perhaps 
because  planktonic  organisms  qualify  as  excellent  candidates  for  being  at  the  mercy  of  their 
environment.  In  a  review  of  the  literature  on  physical  processes  and  planktonic  ecosystems, 
Denman  and  Powell  (1984)  give  numerous  examples  of  successful  results  with  this  approach; 
they  point  out,  however,  that  ecological  responses  often  cannot  be  linked  to  a  particular 
physical  scale.  One  possible  explanation  is  nonlinearity.  Only  in  linear  systems  the  scales  of 
the  response  typically  match  the  scales  of  the  forcings.  Cross-correlation  and  cross-spectral 
analysis  are  examples  of  extensively  used  methods  that  identify  variability  at  similar  scales. 
These  methods  will  therefore  be  most  successful  at  establishing  cause-effect  relationships  in 
linear  systems,  or  close  to  equilibria,  where  nonlinear  systems  are  well  approximated  by  lin¬ 
ear  ones.  There  is,  however,  ample  evidence  for  nonlinearity  and  nonequilibrium  dynamics 
in  population  growth,  ecological  interactions,  and  the  response  of  ecosystems  to  perturba¬ 
tions  (Denman  and  Powell,  1984;  Dwyer  and  Perez,  1983;  Dwyer  et  al. ,  1978;  Ellner  and 
Turchin,  1995;  Turchin  and  Taylor,  1992). 

The  general  theme  of  this  thesis  is  the  interplay  of  scales  in  nonlinear  ecological  systems, 
with  an  emphasis  on  the  mismatch  of  scales  between  biological  variables  and  environmental 
forcings  in  the  plankton.  The  approach  is  theoretical:  I  use  simple  models  to  identify 
conditions  leading  to  such  a  mismatch.  Before  presenting  the  specific  research  chapters, 
I  briefly  clarify  some  terms  that  appear  repeatedly  throughout  this  work:  characteristic 
scales,  external  forcings  and  (non)linearity.  I  also  briefly  review  some  known  ecological 
conditions  for  scale  mismatch  of  environment  and  biology. 
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1.1  On  characteristic  scales  of  variability 


One  definition  of  scale,  commonly  used  in  oceanography,  is  that  of  a  period  or  wavelength 
in  the  temporal  or  spatial  variance  of  a  variable.  Because  variance  generally  occurs  at  more 
than  one  scale,  the  term  characteristic  scale  refers  to  the  dominant  period  or  wavelength. 
A  different  but  related  concept  of  scale  is  the  distance  (or  time)  one  has  to  travel  to  see  a 
significant  change  in  the  quantity  of  interest  (Powell,  1989). 

There  are  many  ways  to  determine  the  characteristic  scale  of  data.  Because  the  power 
spectrum  gives  the  distribution  of  variance  as  a  function  of  frequency  or  wavenumber,  a 
large  peak  in  the  spectrum  indicates  the  prevailing  occurrence  of  variance  at  the  associ¬ 
ated  temporal  or  spatial  scale.  Another  common  measure  of  characteristic  scale  uses  the 
autocorrelation  function.  This  function,  also  known  as  the  correlogram,  gives  the  degree 
of  correlation  between  data  at  different  spatial  or  temporal  lags.  The  lag  at  which  the 
correlogram  first  crosses  zero  is  known  as  the  correlation  length  and  gives  a  measure  of 
characteristic  scale.  The  correlation  length  determines  a  significant  change  in  the  quantity 
of  interest  by  a  significant  decrease  in  the  autocorrelation  function. 

Although  not  presented  here,  there  are  other  measures  of  characteristic  scale.  It  is 
interesting  to  note  that  for  stationary  data,  they  can  all  be  related  to  the  correlogram,  and 
therefore,  to  an  intuitive  interpretation  of  scale.  (See  the  Appendix  for  a  brief  description 
of  other  scale  measures  and  their  respective  relationships  to  the  correlogram). 

1.2  A  sketch  of  a  dynamical  system  with  forcing 

Consider  a  variable  of  interest  Ni,  such  as  the  density  of  a  particular  species  or  the  biomass 
of  a  particular  trophic  level,  whose  changes  in  time  depend  on  its  own  value  and  those  of 
other  variables,  denoted  by  Ni  ( i  =  2,  ...,n).  The  temporal  dynamics  of  Ni  can  be  modelled 
with  a  system  of  (differential)  equations 

fi(NuN2,...,Nn) 
f2(NuN2,....,Nn) 

(1.1) 

fn(NuN2,....,Nn), 
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dNi 

dt 

dN2 

dt 


dNn 

dt 


where  the  functions  /,•  specify  the  respective  rates  of  change  of  the  variables  Nt.  Because 
these  rates  have  no  explicit  dependence  on  time,  system  1.1  has  no  external  forcing  (tech¬ 
nically,  it  is  an  autonomous  system).  Now,  if  Ei(t)  (i  =  l,...,n)  denote  temporal  variables 
affecting  the  rates  /,-,  system  1.1  becomes 


dN  i 
dt 
dN2 

dt 


dNn 

dt 


..,^,^(0) 

=  /2(iVi,JV2,.. 

...,Nn,E2(t)) 

=  fn(NUN2,. 

...,Nn,En(t)) 

(1.2) 


(a  nonautonomous  set  of  equations).  The  external  forcings  Ei(t)  influence  the  dynamics 
without  being  altered  by  the  state  of  the  system.  When  stochastic,  they  are  known  as  dy¬ 
namic  noise;  when  deterministic,  they  can  represent  trends  or  periodic  patterns.  Examples 
of  external  forcings  include:  disturbances  such  as  mortality  due  to  storms,  seasonal  changes 
in  parameters  such  as  temperature  or  mixed  layer  depth,  and  factors  internal  to  the  system 
that  cannot  be  predicted  from  the  state  variables,  such  as  demographic  stochasticity  (Ell- 
ner  and  Turchin,  1995).  The  terms  Ei{t)  are  also  called  exogenous  components,  to  contrast 
them  with  the  endogenous  structure  of  the  system  (equations  1.1)  which  contains  only  the 
feedbacks  between  state  variables  (Ellner  and  Turchin,  1995). 

This  sketch  has  used  differential  equations  and  temporal  dynamics.  The  terminology 
extends,  however,  to  any  type  of  dynamical  formulation  and  to  dynamics  in  space-time. 
Spatial  forcings  are  also  known  as  environmental  heterogeneity. 

Systems  1.1  and  1.2  are  said  to  be  nonlinear  when  the  rate  functions  /,•  are  nonlinear  in 
the  state  variables  Nj.  Most  ecological  models  of  species  interactions  are  nonlinear.  This 
has  important  implications  for  their  response  to  external  forcings. 


1.3  The  mismatch  of  scales  in  ecological  systems:  examples 

Simple  models  can  identify  conditions  that  lead  to  scales  of  variability  in  ecological  patterns 
different  from  those  in  the  underlying  environment.  One  well  known  example,  involving 
nonlinear  ecological  interactions,  is  given  by  predator-prey  models  under  periodic  forcing. 
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These  models  have  revealed  a  rich  array  of  possible  dynamics,  including  frequency-locking, 
quasiperiodicity  and  chaos  (Inoue  and  Kamifukumoto,  1984;  Kot  et  al. ,  1992;  Rinaldi  et  al, 
1993;  Schaffer,  1988).  In  these  dynamic  regimes,  predator  and  prey  can  display  variability 
at  frequencies  other  than  that  of  the  seasonal  forcing.  For  instance,  when  frequency-locked, 
predator  and  prey  solutions  are  periodic  at  multiples  of  the  forcing  period.  A  more  inter¬ 
esting  situation  arises  for  quasiperiodicity:  solutions  are  aperiodic  with  multiple  peaks  in 
the  power  spectrum  at  linear  combinations  of  a  finite  number  of  frequencies,  the  so-called 
fundamental  frequencies  (Parker  and  Chua,  1989).  In  the  predator-prey  models,  quasiperi- 
odic  behavior  involves  two  fundamental  frequencies.  The  scales  associated  with  the  largest 
peaks  in  the  spectrum  may  differ  from  that  of  the  seasonal  forcing.  An  even  more  drastic 
transfer  of  variability  is  exemplified  by  chaos:  solutions  are  aperiodic  with  all  frequencies 
present  in  the  power  spectrum. 

A  critical  condition  for  complex  dynamics  in  these  models  is  the  oscillatory  behavior  of 
predator  and  prey  in  the  absence  of  any  forcing.  These  oscillations  may  be  either  limit  cycles 
or  transient  fluctuations  with  slow  damping.  In  other  words,  the  endogenous  predator-prey 
system  must  have  a  natural  frequency,  and  chaos,  quasiperiodicity  and  frequency-locking 
result  from  the  interplay  of  this  natural  frequency  with  the  seasonal  forcing.  It  is  interesting 
to  note  that  the  propensity  of  planktonic  predator-prey  interactions  to  oscillate  has  been 
observed  both  in  the  field  (McCauley  and  Murdoch,  1987)  and  in  laboratory  experiments 
(Goulden  and  Hornig,  1980;  Pratt,  1983).  Moreover,  complex  dynamics  have  been  found  in 
seasonally  forced  models  for  simple  chemostat  food  webs  (Kot  et  al.,  1992)  and  for  marine 
and  freshwater  planktonic  food  webs  (Caswell  and  Redish,  unpublished;  Doveri  et  al,  1993). 

Predator-prey  models  show  that  a  temporal  scale  mismatch  is  likely  to  occur  when  an 
endogenous  cycle  interacts  with  an  exogenous  frequency.  More  recently,  models  for  epi¬ 
demics  have  revealed  a  more  elaborate  interplay  of  endogenous  cycles  with  periodic  and 
stochastic  forcings  (Engbert  and  Friedhelm,  1994;  Rand  and  Wilson,  1992;  Sidorowich, 
1992).  In  these  nonlinear  models  one  parameter,  the  contact  rate,  varies  seasonally.  The 
stochastic  forcing  results  from  either  low  population  numbers  or  environmental  noise.  In 
the  absence  of  forcing,  the  attractor  of  the  system  is  a  limit  cycle:  the  long-term  solutions 
are  periodic.  With  a  seasonal  contact  rate  but  no  stochasticity,  the  limit  cycle  can  coexist 
in  phase-space  with  a  fascinating  structure  known  as  a  repellor.  Repellors  represent  the  un¬ 
stable  counterparts  of  the  more  familiar  strange  attractors,  that  is,  the  geometrical  objects 
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in  phase-space  onto  which  chaotic  solutions  relax  as  transients  die  out.  Solutions  of  the 
epidemiological  model  are  attracted  towards  the  stable  limit  cycle  but  continuously  pushed 
away,  against  the  unstable  repellor,  by  the  stochastic  forcings  (Rand  and  Wilson,  1992).  In 
this  way,  solutions  continuously  switch  between  short-term  periodic  episodes,  determined 
by  the  limit  cycle,  and  chaotic  transients,  revealing  the  shadow  of  an  unstable  invariant  set. 
These  transients  can  be  long  lasting  because  trajectories  take  a  long  time  to  escape  from 
the  vicinity  of  the  repellor.  In  addition,  because  the  repellor  influences  these  transients, 
solutions  appear  irregular  and  exhibit  variability  at  a  variety  of  scales  not  present  in  the 
environmental  forcings.  This  outcome  may  be  common  in  systems  where  attractors  and 
repellors  interact  with  stochasticity  and  coexist. 


1.4  A  glance  at  what  is  yet  to  come 

The  work  of  the  following  chapters  investigates  with  simple  models  some  novel  scenarios  for 
scale  mismatch  between  environmental  and  ecological  variables.  The  models  are  motivated 
by  planktonic  systems  and  focus  on  one  ubiquitous  nonlinear  ecological  interaction,  that  be¬ 
tween  a  consumer  and  its  resource.  Consumer-resource  interactions  have  been  investigated 
extensively  for  their  response  to  temporal  forcings  (see  section  1.3  for  some  references). 
However,  those  studies  have  for  the  most  part  ignored  the  spatial  dimension  and  the  popu¬ 
lation  structure  generated  by  Me  histories,  two  fundamental  elements  of  ecological  systems 
(Caswell,  1989;  Kareiva,  1994).  Here,  I  extend  the  study  of  consumer-resource  interactions 
under  environmental  forcing  in  two  main  directions:  the  importance  of  population  structure 
and  the  response  to  spatial  heterogeneity.  I  have  organized  this  work  in  three  main  parts 
as  described  below. 

•  Chapters  2  and  3  consider  a  consumer-resource  interaction  with  a  structured  con¬ 
sumer  population.  The  goal  is  to  determine  the  significance  of  population  structure 
to  consumer-resource  dynamics  when  the  resource  supply  varies  in  time.  The  im¬ 
portance  of  population  structure  is  measured  by  its  potential  to  generate  consumer 
patterns  that  differ  in  scale  from  those  of  the  fluctuating  nutrient  supply. 

Specifically,  I  consider  the  interaction  of  a  phytoplankton  population  and  a  limit¬ 
ing  nutrient  resource.  These  interactions  are  often  modelled  with  resource-consumer 
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models  of  the  form 


dN_ 

dt 

dS_ 

dt 


=  YVm 


K  +  S 
=  D(Si-S)- 


N  -DN 
S 


K  +  S 


N 


(1.3) 


These  equations  describe  an  experimental  system  known  as  the  chemostat,  but  have 
been  applied  as  a  basic  building  block  in  planktonic  food  web  models  (Dugdale,  1967; 
Steele  and  Henderson,  1981;  Walz,  1993).  The  variables  N  and  S  denote  phytoplank¬ 
ton  numbers  (biomass  or  density)  and  ambient  resource  concentration,  respectively. 
The  parameters  Vm  and  K  are  the  maximum  rate  and  the  half-saturation  constant  of 
uptake,  Y  is  a  yield  coefficient  converting  units  of  resource  into  units  of  population 
numbers,  D  is  the  dilution  rate  of  the  chemostat,  and  S{  is  the  inflowing  nutrient  con¬ 
centration.  Notice  that  equations  1.3  consider  all  cells  as  equal  and  group  them  under 
a  single  variable  (iV),  the  total  biomass  or  cell  density.  They  do  not  take  into  ac¬ 
count  any  population  heterogeneity  resulting  from  the  life  histories  of  individual  cells. 
However,  single  cells  do  have  life  histories:  each  cell  progresses  through  a  determinate 
sequence  of  events  preceding  cell  division  and  the  population  is  distributed  in  stages 
of  the  cell  cycle.  I  hypothesize  that  the  population  structure  (the  stages  of  the  cell 
cycle)  can  generate  oscillatory  dynamics  in  the  absence  of  a  fluctuating  environment, 
and  that  such  endogenous  oscillations  are  important  to  the  population  response  to 
environmental  variability.  In  Chapter  2, 1  investigate  this  possibility  by  incorporating 
population  structure,  first,  in  a  simple  model  such  as  1.3  ,  and  second,  in  an  extension 
of  the  model  that  allows  cells  to  store  nutrients.  I  explore  the  responses  of  the  models 
to  a  variable  nutrient  supply.  In  Chapter  3,  for  the  purpose  of  comparison,  I  establish 
the  response  of  a  well  known  unstructured  model,  an  extension  of  equations  1.3  with 
nutrient  storage  by  the  cells. 

•  Chapters  4  and  5  add  the  spatial  dimension  to  a  consumer-resource  interaction.  The 
main  goal  is  to  determine  the  consequences  of  environmental  heterogeneity  to  the 
spatio-temporal  dynamics  of  a  predator-prey  interaction.  In  particular,  I  ask:  how 
similar  are  the  population  patterns  to  those  of  the  underlying  environment;  and 
how  does  this  similarity  vary  with  the  type  of  spatio-temporal  dynamics  (periodicity, 
quasiperiodicity  and  chaos). 
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There  are  two  main  classes  of  spatial  models  for  interacting  populations  according 
to  their  treatment  of  space.  Discrete  models  partition  the  environment  into  two  or 
more  patches  whose  dynamics  are  coupled  by  dispersal.  The  second  class  of  models 
represents  space  as  a  continuous  variable.  It  includes  reaction-diffusion  models,  in 
which  time  and  population  densities  are  also  continuous.  Those  models  describe  the 
local  ‘reaction’  of  individuals  (i.e.  the  local  population  dynamics)  and  the  movement 
of  organisms  by  diffusion  (Okubo,  1980).  Reaction- diffusion  equations  are  used  in 
oceanography  to  model  planktonic  systems  in  turbulent  flows  (Okubo,  1980). 

Here,  I  consider  a  reaction-diffusion  equation  for  the  dynamics  of  a  predator  and 
its  prey  in  a  heterogeneous  environment.  I  choose  a  spatial  gradient  as  a  type  of 
environmental  heterogeneity  ubiquitous  in  aquatic  environments  (Mackas  et  al.,  1985). 
In  Chapter  4,  I  focus  on  the  temporal  dynamics  of  the  system;  in  Chapter  5,  on  the 
spatial  consequences  of  the  environmental  gradient. 

•  Chapter  6  stands  alone,  a  little  bit  as  an  outcast,  considering  another  problem  in 
the  approach  of  matching  dominant  scales  of  variability.  Because  nonlinear  systems 
transfer  variability  across  scales,  they  may  lack  a  characteristic  or  dominant  scale.  A 
common  property  of  systems  lacking  a  characteristic  scale  is  scale  similarity:  broadly 
speaking,  a  part  of  the  pattern  resembles  the  whole,  and  therefore,  features  at  one 
scale  are  related  to  those  at  another  by  means  of  one  (or  several)  scaling  factors. 
These  patterns  are  presently  well  known  by  the  name  of  fractals. 

In  spite  of  their  elaborate  forms,  fractals  can  be  described  by  a  single  power  law 
that  relates  some  geometrical  quantity  to  the  scale  at  which  it  is  measured.  They 
are  essentially  geometrical  objects.  Fractals  apply,  however,  to  numerical  data  by 
considering  the  geometrical  properties  of  the  curve  or  surface  associated  with  it  (for 
fractal  descriptions  of  numerical  data,  see  Chapter  4  in  Hastings  and  Sugihara,  1993). 
The  characterization  of  numerical  quantities  that  vary  intermittently  has  motivated 
a  generalization  of  fractals  known  as  multifractals.  Multifractals  describe  patterns 
by  scaling  relations  that  require  a  family  of  different  exponents,  instead  of  the  single 
exponent  of  simple  fractals.  They  have  been  applied  to  a  variety  of  intermittent  mea¬ 
sures  associated  with  nonlinear  phenomena  in  physics  and  geophysics  (see  Sreenivasan, 
1991,  for  a  review  in  fluid  turbulence). 
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In  Chapter  6, 1  give  a  first  application  to  biological  oceanography.  I  analyze  an  acous¬ 
tic  data  set  on  zooplankton  biomass  to  describe  the  distribution  in  time  of  the  total 
variability  in  the  data.  This  distribution  is  highly  intermittent:  extreme  localized  con¬ 
tributions  account  for  a  large  proportion  of  total  variability.  I  show  that  multifractals 
provide  a  good  characterization  of  such  variability. 

After  the  research  chapters,  a  short  section  recapitulates  some  main  results  in  light  of 
the  general  theme  presented  here. 

Biological  oceanography  has  pioneered  the  concept  of  scale  in  ecology  (Haury  et  al. ,  1978; 
Steele,  1978),  but  a  linear  perspective  has  dominated  the  view  of  how  environmental  and 
biological  scales  interact.  A  few  authors  have  cautioned  against  simple  linearity  assumptions 
(Denman  and  Powell,  1984;  Dwyer  and  Perez,  1983;  Star  and  Cullen,  1981;  Steele,  1988). 
The  following  chapters  should  further  support  the  need  for  a  nonlinear  perspective. 
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Chapter  2 


Phytoplankton  population 
dynamics  and  nutrient  variability: 
a  cell  cycle  perspective 


...each  cell  is  a  citizen... 

— Schwann,  trans.  H.  Smith  in  Schwann  and  Schleiden 
Researches ,  1847. 


2.1  Introduction 

The  dynamics  of  nonlinear  systems  depend  on  the  interplay  of  time  scales.  In  population 
interactions,  one  of  these  time  scales  is  provided  by  the  process  of  reproduction,  growth  and 
maturation;  i.e.,  by  the  generation  time  of  the  organism.  Structured  population  models 
classify  individuals  by  age,  size,  or  developmental  stage  to  incorporate  this  time  scale. 
Unstructured  models,  written  in  terms  of  bulk  variables  such  as  total  numbers  or  biomass 
do  not. 

Even  unicellular  organisms  like  phytoplankton  have  a  life  history:  the  cell  division  cycle. 
Each  cell  progresses  through  a  cyclic  sequence  of  events  preceding  cell  division  (Mitchison, 
1971),  and  the  population  is  distributed  among  different  stages  of  the  cell  cycle  (Figure  2-1). 
This  population  heterogeneity  is  potentially  important  to  population  dynamics  because  it 
interacts  with  the  environment:  resource  levels  may  affect  both  uptake  by  the  cells  and 
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cell  progression  through  the  cycle  at  specific  stages  (Brezinski,  1992;  Olson  and  Chisholm, 
1986;  Vaulot  et  al.,  1987).  Experiments  have  shown  that  in  different  phytoplankton  species, 
nutrient  deprivation  blocks  progression  through  the  cycle  at  specific  stages,  and  that  the 
duration  of  such  stages  also  lengthens  under  nutrient  limitation  (Olson  and  Chisholm,  1986, 
Olson  et  al.  1986;  Vaulot,  1985;  Vaulot  et  al.,  1987).  These  results  are  consistent  with  a 
conceptual  view  of  the  cell  cycle,  known  as  the  transition  point  hypothesis,  in  which  an 
environmental  factor  has  no  effect  on  cell  progression  beyond  a  certain  point  in  the  cycle 
(Spudich  and  Sager,  1980;  Vaulot  et  al.,  1986). 

Two  main  environmental  factors,  light  and  nutrients,  control  progression  of  a  phyto¬ 
plankton  cell  through  its  cycle  (Prezelin,  1992).  Population  models  have  shown  that  the 
light-dark  cycle  coupled  to  the  transition  point  hypothesis  generates  and  explains  observed 
oscillations  in  division  patterns  (Heath  and  Spencer,  1985;  Vaulot,  1985).  While  the  con¬ 
sequences  of  the  photoperiod  on  population  dynamics  through  the  cell  cycle  are  therefore 
well  understood,  those  of  nutrient  fluctuations  are  not.  Because  nutrients  are  consumed  by 
the  cells,  there  is  however  a  fundamental  difference  between  these  two  driving  forces  of  the 
cell  cycle. 

To  investigate  the  significance  of  the  cell  cycle  to  the  dynamics  of  phytoplankton  pop¬ 
ulations,  this  chapter  describes  a  chemostat  model  that  incorporates  the  population  distri¬ 
bution  along  the  cell  cycle  and  the  transition  point  hypothesis.  I  will  show  that  under  a 
constant  supply  of  nutrients  (and  light),  total  cell  numbers  are  capable  of  oscillatory  dynam¬ 
ics.  These  oscillations  are  generated  by  the  interplay  between  the  environmental  resource 
levels  and  the  population  distribution  in  stages  of  the  cell  cycle.  They  introduce  a  bio¬ 
logical  frequency  capable  of  interacting  with  environmental  forcing  frequencies  to  generate 
complex  temporal  dynamics.  Under  a  periodic  nutrient  supply,  total  cell  numbers  display 
aperiodic  dynamics  with  variability  at  frequencies  other  than  that  of  the  forcing.  I  then 
extend  the  model  to  incorporate  nutrient  storage  by  the  cells.  Cell  progression  through  part 
of  the  cycle  becomes  a  function  of  the  cellular  nutrient  levels.  The  main  qualitative  results 
on  population  dynamics  remain  unchanged.  I  discuss  how  they  differ,  however,  from  the 
dynamics  of  traditional  chemostat  models  and  from  cell  cycle  models  driven  by  the  light- 
dark  cycle.  Finally,  I  propose  a  link  between  resource  control  of  cell  cycle  progression  and 
the  time  delays  between  resource  levels  and  population  growth,  previously  postulated  to 
explain  oscillatory  transients  in  chemostat  experiments  (Caperon,  1969;  Cunningham  and 
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Figure  2-1:  The  cell  cycle.  The  cell  cycle  is  classically  divided  into  four  stages.  The  cell 
genome  replicates  during  S,  M  corresponds  to  mitosis  and  cell  division,  G1  and  G2  denote 
stages  during  which  most  of  cell  growth  takes  place.  Adapted  from  Chang,  1989. 

Maas,  1978;  Cunningham  and  Nisbet,  1980;  Williams,  1971). 

Cell  cycle  modelling  for  non-marine  cells  is  an  active  field  of  research,  but  has  developed 
with  an  emphasis  on  cellular  as  opposed  to  ecological  phenomena.  In  the  last  decade,  phy¬ 
toplankton  ecology  has  moved  into  the  small  scales  of  the  individual  both  in  the  laboratory 
and  in  the  field  (Chisholm  et  al.,  1986;  Heath,  1988;  Harris,  1980).  It  is  an  open  question 
whether  these  small  scales  influence  the  dynamics  at  higher  levels. 

2.2  The  basic  model 

The  setting  for  the  model  is  the  method  of  continuous  culture  known  as  the  chemostat.  This 
choice  allows  comparisons  of  the  model  dynamics  to  experimental  results  in  the  literature, 
and  to  the  well-known  behavior  of  unstructured  chemostat  models  grouping  all  cells  into  a 
single  variable  (Monod,  1942;  Droop,  1974;  Lange  and  Oyarzun,  1992;  Smith  and  Waltman, 
1994). 

The  chemostat  provides  a  simple,  yet  controllable  idealization  of  an  aquatic  system 
with  both  an  inflow  and  an  outflow  of  nutrients.  Nutrients  at  an  input  concentration  Si, 


are  supplied  by  a  through  flow  at  rate  F  into  a  chamber  of  volume  V .  The  effluent  contains 
both  medium  and  phytoplankton  cells,  and  the  residence  time  of  the  cells  in  the  chamber 
is  given  by  the  reciprocal  of  the  dilution  rate  D  =  F/V. 

2.2.1  Equations 

Equations  for  the  dynamics  of  cell  populations  have  been  formulated  in  both  discrete  and 
continuous  time.  Discrete  models  divide  the  cell  cycle  into  discrete  stages  such  as  the  four 
conventional  stages  G1-S-G2-M  (Figure  2-1),  or  the  two  parts  of  the  cycle  separated  by 
a  transition  point  (Smith  and  Martin,  1973;  Heath  and  Spencer,  1985).  I  choose  here  the 
continuous  representation.  A  variable,  denoted  by  p,  measures  the  extent  of  cell  development 
or  position  along  the  cell  cycle  (Hoppensteadt,  1986;  Eubinow,  1968).  The  rate  of  change 
of  p  with  time,  the  maturation  velocity  v  =  dp/dt ,  describes  cell  progression  through  the 
cycle. 


Figure  2-2:  Conceptual  representation  of  the  cell  cycle.  In  the  model,  a  variable  p  indicates 
the  position  of  a  cell  along  the  cell  cycle.  The  interval  [po,pc]  is  a  nutrient-dependent 
segment  during  which  progression  of  a  cell  through  the  cycle  depends  on  nutrient  levels. 
In  the  rest  of  the  cycle,  cells  progress  at  a  constant  rate.  Thus,  once  a  cell  reaches  pc,  it 
proceeds  towards  division  regardless  of  nutrient  conditions. 

A  variable  N(p,t )  describes  the  distribution  of  cells  along  the  cell  cycle.  Notice  that 
total  cell  numbers  Ntot(t )  is  obtained  by  integrating  this  distribution  over  p, 
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(2.1) 


Nt <>t(t)  =  f  N(p,t)dp. 

JP 

The  dynamics  of  N(p,  t )  is  given  by  the  following  equation 

^  =  —m(p)N  (2.2) 

(an  extension  due  to  Rubinow  (1968)  of  the  Me  Kendrick- Von  Foerster  equation  for  age 
structured  populations  (Von  Foerster,  1959)),  where  m  denotes  the  mortality  rate. 

The  variable  p  is  normalized  so  that  cells  are  born  with  p  —  0  and  the  average  cell 
divides  at  p  =  1.  Although  p  is  continuous,  the  model  incorporates  discrete  stages  by 
subdividing  the  interval  [0, 1]  into  subintervals,  and  specifying  equation  2.2  for  each  of 
these  subintervals.  To  incorporate  the  transition  point  hypothesis,  I  divide  the  cell  cycle 
into  two  different  portions  and  introduce  a  transition  point  pc  between  resource-dependent 
and  resource-independent  segments  (Figure  2-2).  In  the  subinterval  [po,pc]  progression  of 
a  cell  through  the  cycle  is  a  function  of  ambient  nutrient  levels.  More  specifically,  the 
maturation  velocity  in  this  part  of  the  cycle  is  proportional  to  nutrient  uptake, 


dp  S 
~dt  =  V°  K  +  S 


(2.3) 


where  nutrient  uptake  follows  a  Monod  type  curve  with  half-saturation  constant  K,  and  vq 
is  the  maximum  maturation  rate.  Through  the  rest  of  the  cycle,  cells  progress  at  a  constant 


rate  given  by 


(2.4) 


Once  a  cell  reaches  the  transition  point  pc,  it  proceeds  towards  division  regardless  of 
environmental  conditions.  By  specifying  equation  2.2  in  each  part  of  the  cycle  and  using 
equations  2.3  and  2.4,  the  population  model  becomes 


dN  S  ON  .  ,  ,  , 

=  -DN  fotpe[po'pJ  (2'6) 

dN  ON 

— — b  is c ~~p.  =  —DN  -  B(p)N  otherwise  (2-6) 

at  op 


where  the  loss  rate  m  (equation  2.2)  is  replaced  by  the  dilution  rate  D  and  the  division 
rate  B(p).  The  population  model  is  completed  by  describing  cell  division  with  a  boundary 
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condition  for  the  flux  of  newborn  cells  at  p  =  0, 

N( 0,  t)uc  =  2  f  B(p)N(p,  t)dp,  (2.7) 

JP 

were  each  cell  divides  into  two  daughter  cells.  To  specify  the  division  rate  B(p),  I  consider 
that  there  is  population  stochasticity  in  the  value  of  p  at  which  cells  divide,  and  describe 
this  variability  by  4>(p ),  a  probability  density  for  the  maturity  stage  at  division  in  a  cohort 
of  newborn  cells  (i.e.  <j>(p)dp  is  the  proportion  of  cells  dividing  between  maturity  stage  p 
and  p  +  dp).  Once  <f>(p)  is  specified,  the  division  rate  is  obtained  as 

_  dP  4jp) 

dt  (1  —  Jq  <j>(s)ds) 

(Metz  and  Diekman,  1980). 

Finally,  the  dynamics  of  the  ambient  nutrient  concentration  are  given  by 

(2.9) 

where  Vm  denotes  the  maximum  uptake  rate  of  a  cell  and  Si,  the  inflowing  nutrient  con¬ 
centration.  The  ambient  nutrient  concentration  increases  with  the  inflow  of  nutrients  to 
the  chemostat,  and  decreases  with  the  ouflow  of  nutrients  and  with  uptake.  Differences  in 
uptake  among  cells  are  considered  negligible  and  total  uptake  is  computed  by  multiplying 
cell  uptake  by  total  cell  numbers  (see  section  2.5  for  a  discussion  of  this  assumption). 

2.2.2  Numerical  methods  and  non-dimensional  equations 

The  dynamics  of  the  system  is  investigated  by  simulating  the  model  with  the  numerical 
method  known  as  the  Escalator  box-car  train  (De  Roos  et  al.  1992;  De  Roos,  1988).  This 
method  was  developed  for  the  integration  of  partial  differential  equations  modelling  pop¬ 
ulations  structured  by  variables  other  than  age.  It  subdivides  the  population  into  cohorts 
whose  evolution  is  followed,  making  it  possible  to  track  cells  that  have  experienced  the 
same  environment.  For  a  description  of  the  application  of  the  method  to  this  model,  see 
the  Appendix. 

To  reduce  the  number  of  parameters  and  focus  on  the  qualitative  dynamics  of  the  system, 
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I  simulate  the  model  in  non-dimensional  form.  With  the  non-dimensional  variables, 

S  NVm 

s  =  — ,  n=-p — »  T  =  tu o,  (2.10) 

A  Iu/q 

the  population  model  becomes, 

dn  s  dn 

+  =  ~dn  {orpe  [p°'Pc]  (2'n) 

dn  dn  ,  ,  , 

— — f-  u—  =  —  (d  +  b)n  otherwise  (2.12) 

dr  dp 

with  boundary  condition 

tm(0,r)  =  2  f  b(p)n(p,  r)dp,  (2.13) 

JP 

and  resource  dynamics  given  by 

ds  ,  s  .  . 

—  =  d{Si  —  s)  —  ^  ^  ntot ■  (2-14) 

The  non-dimensional  parameters  in  the  above  equations  are  the  dilution  rate  d  =  D/v o,  the 
division  rate  b(p )  =  B(p)/vo ,  the  maturation  velocity  in  the  nutrient-independent  segment 
v  =  ucjv o,  and  the  inflowing  nutrient  concentration  S{  =  Si/K. 

In  the  simulations,  the  maturation  velocities  uo  and  vc  are  assumed  equal  (i.e.  v  =  1). 
The  normal  distribution  is  used  for  the  probability  density  of  the  maturity  stage  at  division 
4>(p).  This  is  consistent  with  observed  bell-shaped  distributions  of  cell  cycle  duration  for  a 
variety  of  microorganisms  under  constant  environmental  conditions  (Cook  and  Cook,  1962; 
Miyata  e  al.,  1978;  Prescott,  1959).  The  mean  of  the  distribution  was  located  at  p  =  1,  and 
the  variance,  a2  =  0.01,  was  selected  sufficiently  small  to  make  division  negligible  before 
the  transition  point  pc.  Two  of  the  parameters,  the  non-dimensional  dilution  rate  d  and 
the  inflowing  nutrient  concentration  s;  are  under  experimental  control.  The  value  of  S{ 
compares  the  inflowing  nutrient  concentration  to  the  half-saturation  constant  of  nutrient 
uptake.  The  dilution  rate  d  compares  the  rate  of  cell  loss  from  the  chemostat  to  the  rate  of 
cell  progression  through  the  cycle.  Equivalently,  it  compares  the  residence  time  of  the  cells 
in  the  chemostat  to  the  mean  time  it  takes  for  a  cell  to  complete  the  cycle  when  resource 
levels  are  high.  Thus,  the  lower  the  value  of  d,  the  longer  cells  stay  in  the  chemostat  chamber 
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with  respect  to  their  cycling  time. 

All  the  figures  on  the  model  dynamics  are  plotted  in  non-dimensional  units.  Parameter 
values  from  the  literature  on  chemostat  experiments  can  be  used  to  convert  the  axes  to 
dimensional  units.  For  instance,  non-dimensional  cell  numbers  multiplied  by  (Ku o)/Vm 
(~  106  to  10s,  DiToro,  1980),  give  cells  per  liter.  Non-dimensional  time  divided  by  uq  (~  1 
to  4  if  literature  values  for  maximum  division  rates  are  used,  Di  Toro  ,  1980)  gives  time  in 
days.  These  ranges  encompass  data  for  a  variety  of  species  and  limiting  nutrients. 

2.2.3  Model  dynamics  under  a  constant  environmental  forcing 

I  consider  first  a  nutrient  supply  that  is  constant  in  time,  and  focus  on  qualitative  changes  in 
dynamics  resulting  from  different  values  of  d  and  s;,  the  two  parameters  under  experimental 
control. 


Time(T) 

Figure  2-3:  Oscillatory  transients  for  increasing  nutrient  inflow.  Total  cell  numbers  converge 
to  a  steady-state  for  s,-  —  1.  With  this  steady-state  as  initial  condition,  the  model  was  run  for 
increasing  values  of  s,-.  Transient  oscillations  appear  (s,-  =  2,3).  For  s;  =  5  the  oscillations 
persist,  (v  =  1,  [po,Pc]  —  [0.1, 0.5],  d  =  0.3.) 
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Simulations  indicate  that  changes  in  these  parameters  modify  the  model  behavior  from 
steady-state  to  oscillatory  dynamics.  This  is  illustrated  in  Figure  2-3  .  When  s;  =  1,  total 
cell  numbers  converge  to  an  equilibrium.  At  equilibrium,  the  population  reaches  a  stable 
state  distribution  (i.e.  a  distribution  along  the  cell  cycle  that  does  not  change  in  time) 
(Figure  2-4  ).  In  Figure  2-3  all  simulations  have  the  same  initial  condition  given  by  the 
equilibrium  for  S{  =  1.  The  model  is  run  for  increasing  values  of  s;.  Oscillatory  transients 
of  increasing  amplitude  appear  (s;  =  2  and  Si  =  3).  For  s;  =  5,  these  oscillations  persist 
and  converge  to  a  limit  cycle  (Figure  2-3  and  2-5(A)).  Similarly,  as  d  decreases  and  the 
residence  time  of  the  cells  in  the  chemostat  increases,  oscillatory  transients  and  persistent 
oscillations  appear.  Figure  2-6  shows  the  behavior  of  total  cell  numbers  for  different  values 
of  d.  Initial  conditions  are  given  by  the  above  limit  cycle  (i.e  s,-  =  5  and  d  -  0.3).  For 
large  d,  the  outflow  of  cells  from  the  chemostat  is  too  fast  for  the  population  to  persist.  For 
d  =  0.5,  oscillations  are  damped  and  the  population  converges  to  a  steady-state.  For  low  d 
values,  limit  cycles  occur  (  Figure  2-5(B)).  As  d  decreases  both  the  amplitude  and  period 
of  the  cycles  increase  (Figure  2-6). 


Figure  2-4:  Normalized  population  distribution.  At  steady-state  the  population  reaches  a 
stable  distribution  along  the  cell  cycle.  The  integral  under  the  curve  gives  the  fraction  of 
the  total  population  in  an  interval  of  p.  The  squares  correspond  to  the  different  cohorts  in 
the  simulation.  They  indicate  the  fraction  of  the  total  population  in  a  cohort  with  mean 
maturity  stage  p.  (s;  =  1,  v  =  1,  [po,Pc]  =  [0.1, 0.5],  d  =  0.3. 
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A 
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Figure  2-5:  Limit  cycles.  Total  cell  numbers  vs.  ambient  nutrient  levels  after  transients 
have  died  out.  In  (A),  d  =  0.3;  in  (B),  d  =  0.1.  (s;  =  5,  v  =  1,  [po>Pc]  =  [0.1, 0.5]). 


Time(T) 


Figure  2-6:  Oscillatory  dynamics  for  low  dilution  rates.  At  low  dilution  rates  the  model 
converges  to  a  limit  cycle  (d  =  0.3  and  d  =  0.1).  For  higher  d  =  0.5,  the  population 
reaches  a  steady-state.  For  d  —  0.7,  cell  losses  are  too  high  for  the  population  to  persist 
in  the  chemostat.  Initial  conditions  are  the  limit  cycle  for  d  =  0.3.  (s;  =  5,  v  =  1, 
[pchPc]  =  [0.1,  0.5]). 
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The  oscillatory  behavior  of  the  model  results  from  the  interaction  between  ambient 
resource  levels  and  the  population  distribution  along  the  cycle.  Figure  2-7  shows  a  complete 
cycle  of  the  ambient  nutrient  level  6,  total  cell  numbers,  and  the  number  of  cells  in  the 
nutrient-dependent  segment  [p0,pc].  As  total  cell  numbers  increase,  ambient  nutrient  levels 
decrease.  Then,  more  and  more  cells  accumulate  in  [po,pc]  as  their  progression  through  this 
part  of  the  cycle  slows  down.  As  cells  accumulate  in  [ Po,Pc ],  fewer  cells  are  able  to  complete 
the  cycle  and  divide,  and,  total  cell  numbers  and  nutrient  uptake  decrease.  Ambient  nutrient 
levels  then  go  up.  This  permits  the  cells  in  [poiPc]  to  proceed  towards  division  and  results 
in  the  next  pulse  in  cell  numbers.  The  whole  cycle  restarts  again.  These  fluctuations,  called 
here  generation  cycles,  differ  from  typical  predator-prey  cycles  in  which  the  whole  predator 
population  would  oscillate  in  synchrony  without  changes  in  population  structure. 


Figure  2-7:  The  anatomy  of  generation  cycles.  A  complete  oscillation  of  the  model  is  shown 
for  total  cell  numbers,  cell  numbers  in  the  segment  [po,pc],  and  ambient  nutrient  levels.  The 
model  oscillations  involve  changes  in  the  population  distribution.  For  a  complete  description 
of  this  figure  see  text  (section  2.2.3).  (st  =  5,  v  =  1,  [po,Pc]  =  [0.1, 0.5],  d  =  0.3). 
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2.2.4  Model  dynamics  under  a  variable  nutrient  supply 

The  generation  cycles  introduce  a  characteristic  population  frequency  capable  of  interact¬ 
ing  with  external  environmental  frequencies  to  generate  complex  temporal  dynamics.  To 
investigate  this  possibility,  I  consider  the  simplest  form  of  a  variable  nutrient  supply,  the 
periodic  function 


£(1  +  esin  (— <)) 


(2.15) 


whith  mean  value  5,,  amplitude  e  and  period  T. 

With  this  periodic  input,  the  non-dimensional  equation  for  5  =  S/K  becomes 


ds 

dr 


d(si(  1  +  csinwr)  —  s ) 


s 

1  +  3 


^tot' 


(2.16) 


where  the  frequency  u>  —  (2tt)/(voT)  and  st  =  Si/K.  Other  equations  remain  the  same. 

The  model  exhibits  two  types  of  response  to  a  periodic  nutrient  supply:  periodic  and 
aperiodic  dynamics.  It  is  the  latter  that  I  wish  to  emphasize  here  since  this  aperiodic 
dynamics  exhibits  variability  at  frequencies  other  than  that  of  the  forcing.  Thus,  the  popu¬ 
lation  is  capable  of  a  more  complex  response  than  a  simple  cycle  tracking  the  environmental 
forcing.  Figures  2-8  and  2-9  illustrate  the  model  dynamics  for  two  different  forcing  frequen¬ 
cies,  oj  =  4.5  and  u  =  3.2  respectively.  For  comparison,  the  natural  frequency  of  the  system 
(i.e.  the  frequency  of  the  generation  cycles)  is  w  =  2.5.  The  behavior  of  total  cell  numbers 
after  transients  have  died  out  is  aperiodic  (Figures  2-8(A)  and  2-9(A)).  This  can  be  seen 
by  plotting  one  of  the  variables,  for  instance  the  nutrient  concentration  s(r),  vs.  itself  at 
lagge<i  intervals  of  time  (Figure  2-10).  If  the  dynamics  were  periodic  the  trajectory  would 
come  back  on  itself.  Instead,  the  trajectory  moves  on  the  surface  of  a  torus  and  never 
repeats  itself.  This  behavior  is  known  as  quasiperiodic  dynamics.  By  contrast  to  periodic 
oscillations,  which  have  only  one  fundamental  frequency,  quasiperiodic  behavior  has  two  or 
more  fundamental  frequencies  (two  in  the  case  of  a  torus  attractor).  Its  power  spectrum  can 
display  peaks  at  harmonics  of  these  fundamental  frequencies  and  at  sums  and  differences 
of  these  harmonics.  Figure  2-11  shows  the  power  spectrum  of  the  solution  for  total  cell 
numbers.  The  arrows  indicate  the  two  dominant  frequencies.  In  this  particular  example, 
one  of  these  frequencies  coincides  with  that  of  the  forcing,  the  other  one,  to  that  of  the 
generation  cycles.  This  is  not  always  the  case.  As  the  frequency  of  the  forcing  approaches 
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that  of  the  generation  cycles,  one  of  the  two  dominant  frequencies  can  differ  from  both.  For 
instance,  the  low  frequency  modulation  of  the  solution  in  Figure  2-9  is  not  present  in  the 
nutrient  forcing,  and  is  lower  than  the  natural  beating  of  the  system. 
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Figure  2-8:  Quasiperiodic  dynamics  under  a  periodic  nutrient  supply.  Total  cell  numbers 
(A)  exhibit  an  aperiodic  response  to  the  periodic  nutrient  supply  (B).  Only  the  long  term 
population  behavior  is  shown,  (s,  =  5,  v  =  1,  [po,pc]  =  [0.1, 0.5],  d  =  0.3,  e  =  0.9,  u  =  4.5). 
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Figure  2-9:  Quasiperiodic  dynamics  under  a  periodic  nutrient  supply.  The  periodic  nutrient 
inflow  is  shown  in  (B),  the  aperiodic  response  of  the  population,  in  (A).  (st-  =  5,  v  =  1, 
[po,  pc]  =  [0.1, 0.5],  d  =  0.3,  e  =  0.9,  u  =  3.2). 
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Figure  2-10:  Torus  attractor.  The  attractor  of  the  system  is  reconstructed  by  plotting 
one  of  the  variables,  the  ambient  nutrient  s,  against  itself  at  lagged  intervals  of  time  after 
transients  have  died  out.  The  trajectory  moves  on  the  surface  of  a  torus,  (s,-  =  5,  v  =  1, 
[Po,Pc]  =  [0.1, 0.5],  d  =  0.3,  e  =  0.9,  w  =  4.5). 


Frequency  (27ix)’ 

Figure  2-11:  Power  spectrum  of  quasiperiodic  behavior.  The  spectrum  of  total  cell  numbers 
in  Figure  2-8  shows  variance  at  dominant  frequencies  u  (the  forcing  frequency)  and  ugc  (the 
natural  frequency  of  the  generation  cycles).  ( s;  =  5,  v  =  1,  \p0,pc ]  =  [0.1, 0.5],  d  =  0.3, 
<r  =  0.9,  u  =  4.5). 
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Figure  2-12:  Irregular  transients  in  total  cell  numbers.  For  a  periodic  nutrient  supply 
total  cell  numbers  display  irregular  fluctuations.  In  (A),  the  transients  correspond  to  the 
simulation  shown  in  Figure  2-9;  in  (B),  to  the  simulation  shown  in  Figure  2-8. 

The  above  results  pertain  to  the  long  term  behavior  of  the  system.  Transient  dynamics, 
which  are  more  relevant  to  chemostat  experiments,  take  a  long  time  to  die  out.  They  share, 
however,  many  of  the  properties  of  the  long  term  dynamics  (Figure  2-12).  They  are  aperi¬ 
odic  and  do  not  simply  track  the  environmental  forcing.  As  a  result,  the  cross- correlation 
between  population  numbers  and  nutrient  forcing  is  low  for  any  time  lag  (Figure  2-13). 
Thus  observations  of  such  a  system  would  suggest  only  a  weak  link  between  phytoplankton 
and  nutrient  input. 

For  some  forcing  frequencies,  the  response  of  the  model  is  periodic.  An  example  is  shown 
in  Figure  2-14  .  Although  total  cell  numbers  oscillate  at  the  environmental  frequency,  they 
display  multiple  peaks  within  a  cycle  which  are  not  present  in  the  forcing. 
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Figure  2-14:  Periodic  response  to  periodic  forcing.  Total  cell  numbers  oscillates  with  the 
forcing  frequency  but  displays  multiple  peaks  in  each  cycle. 
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2.3  Population  model  with  nutrient  storage  by  the  cells 


Phytoplankton  cells  are  known  to  be  capable  of  storing  nutrients.  Because  nutrient  storage 
temporally  decouples  cell  physiology  from  ambient  nutrient  levels,  it  plays  an  important 
role  in  variable  environments  (Harris,  1980).  Chemostat  models,  introduced  as  extensions 
of  the  classical  Monod  equations  (Monod,  1942),  incorporated  nutrient  storage  by  the  cells 
and  made  population  growth  a  function  of  internal  nutrient  levels  (Droop,  1973;  Caperon, 
1969).  Here,  I  extend  the  above  cell  cycle  model  to  consider  an  intracellular  store  of  nutrient. 
Progression  through  the  nutrient-dependent  segment  of  the  cell  cycle  becomes  a  function 
of  internal  nutrient  levels.  With  this  model  extension,  I  ask  if  nutrient  storage  by  the  cells 
modifies  the  qualitative  results  obtained  with  the  original  model. 


2.3.1  Equations 

An  additional  variable,  the  cell  quota  Q,  is  introduced,  and  defined  as  the  amount  of  stored 
nutrient  in  a  cell.  The  population  density  N(t,p,Q)  becomes  a  function  of  both  p  and  Q. 

I  specify,  first,  the  dynamics  of  the  population  in  the  nutrient-dependent  segment  [po,pc]- 
The  dynamics  of  the  population  density  N(t,  p,  Q)  is  given  by  the  following  balance  equation 


dN  dpdN  dQ  dN 
dt  ^  dt  dp  ^  dt  dQ 


(2.17) 


in  which  the  total  derivative  of  N  equals  cell  losses  due  to  outflow.  Through  this  part  of 
the  cycle,  the  cell  quota  determines  the  maturation  speed  of  a  cell.  Thus, 


dp 

dt 


«/o(l-Tf)  if  Q  >  Kq 
0  otherwise 


(2.18) 


where  Kq  is  a  treshold  below  which  cell  progression  through  the  cycle  stops.  The  expres¬ 
sion  for  the  maturation  speed  in  equation  2.18  is  borrowed  from  the  function  proposed  by 
Droop  (1974)  for  population  growth.  It  has  two  main  properties:  saturation  at  a  maximum 
value  of  vq  and  the  existence  of  a  treshold  for  growth. 

In  the  rest  of  the  cycle,  cells  proceed  at  a  constant  speed  vc.  The  population  equation 
is  then  given  by 


dN  dN  dQ  dN 
dt  ^  Vc  dp  dt  dQ 


~(B(p)  +  D)N 


(2.19) 
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where  cell  losses  occur  both  by  division  and  outflow  from  the  chemostat. 

The  cell  quota  increases  due  to  nutrient  uptake  and  decreases  by  cell  division,  which 
partitions  the  stored  nutrients  equally  bewteen  two  daughter  cells.  Then,  the  rate  of  change 
of  Q  for  any  cell  that  is  not  dividing  is  given  by 


dQ 

dt 


(2.20) 


where  Vm  and  Ku  denote,  respectively,  the  maximum  uptake  rate  and  the  half-saturation 
constant  of  uptake.  As  before,  nutrient  uptake  is  a  saturating  function  of  the  ambient 
nutrient  concentration. 

Cell  division  is  specified  as  a  boundary  condition  giving  the  flux  of  cells  at  p  =  0 


vcN(t,0,Q/2)  =  2  [  B(p)N(t,p,Q)dp.  (2.21) 

JP 

The  model  is  completed  with  the  equation  for  the  dynamics  of  the  ambient  nutrient  con¬ 
centration, 

ft  (2.22) 

As  before,  uptake  differences  among  cells  are  considered  negligible  (see  section  2.5  for  a 
discussion  of  this  assumption). 


2.3.2  Numerical  method  and  non-dimensional  equations 

Simulations  with  the  Escalator  boxcar  train  method  follow  the  number  of  cells  in  different 
cohorts,  as  well  as  the  mean  maturation  stage  and  the  mean  cell  quota  for  each  cohort  (see 
Appendix).  To  reduce  the  number  of  parameters,  the  model  was  written  in  nondimensional 
form.  The  non-dimensional  variables  are 


T  =  fl/Q. 


Then,  the  population  dynamics  is  given  by 


dn  dp  dn 
dr  +  dr  dp 


s  dn 
1  +  s  dq 


=  —dn 


(2.23) 


(2.24) 


40 


for  p  in  [po,Pc],  where 


dp 

dr 


(1-i)  if?>l 

<  v 

0  otherwise. 


(2.25) 


In  the  rest  of  the  cycle, 


dn 

dr 


dn 

+  v—  +  U 
op 


s  dn 
l  +  s~dq 


—(d  +  b)n 


(2.26) 


with  boundary  condition 


vn(T,0,q/2)  =  2  f  b(p)n(r,p,q)dp.  (2.27) 

Jp 

The  resource  dynamics  are  given  by 

^.  =  d(si-s)-U-^-ntot.  (2-28) 

dr  1  +  5 

The  non-dimensional  parameters  are  the  dilution  rate  d  =  D/oq,  the  division  rate  b(p)  = 
B(p)/v o,  the  maturation  velocity  in  the  nutrient-independent  segment  v  =  vc/v o,  the  in¬ 
flowing  nutrient  concentration  s;  =  S\/Ku,  and  the  maximum  uptake  rate  U  =  Vm/(KQU0). 
The  parameter  U  characterizes  a  phytoplankton  species  with  respect  to  a  particular  nutri¬ 
ent  by  comparing  the  temporal  scales  of  population  growth  and  nutrient  uptake.  If  vq  is 
estimated  from  literature  values  on  maximum  division  rates,  experimental  work  for  various 
phytoplankton  species  shows  that  U  spans  a  broad  range  (U  ~  50  -  200  for  phosphorus, 
U  ~  5  -  25  for  nitrogen,  U  ~  1  -  3  for  silica)  (DiToro,  1980).  To  convert  the  figures  axes 
to  dimensional  units  divide  time  by  u0  (~  1  -  4  ,  if  literature  values  for  maximum  division 
rates  are  used)  to  obtain  time  in  days,  and  multiply  cell  numbers  by  (~  10'  -  109  for 
different  phytoplankton  species  and  limiting  nutrients  (DiToro,  1980)),  to  obtain  cells  per 
liter. 

2.3.3  Model  dynamics  under  a  constant  nutrient  supply 

I  focus  again  on  the  qualitative  changes  in  the  model  dynamics  as  the  two  parameters  under 
experimental  control,  d  and  st,  are  varied. 

The  main  results  for  the  qualitative  dynamics  of  the  original  model  hold  when  nutrient 
storage  by  the  cells  is  considered.  The  population  is  capable  of  oscillatory  dynamics.  These 
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cycles  occur  for  low  values  of  the  dilution  rate,  that  is,  for  long  residence  times  of  the  cells 
in  the  system.  For  high  dilution  rates,  the  population  reaches  an  equilibrium  and  a  stable 
distribution  (Figures  2-15  and  2-16,  d  =  0.5).  As  the  dilution  rate  decreases,  oscillatory 
transients  appear  and  eventually  the  equilibrium  becomes  unstable  and  the  system  converges 
to  a  limit  cycle.  Figure  2-16  illustrates  the  changes  in  dynamics  for  different  dilution  rates. 
For  d  <  0.2,  transient  oscillations  of  small  amplitude  decay  to  steady-state.  For  d  =  0.1 
the  oscillations  persist  and  converge  to  a  limit  cycle  (Figure  2-17  ).  Note  also  the  long 
transients  characteristic  of  this  model. 


Figure  2-15:  Normalized  population  distribution.  At  steady-state  the  population  reaches  a 
stable  distribution.  The  integral  under  the  curves  gives  the  fraction  of  thedotal  population 
in  an  interval  of  p  (A)  or  q  values  (B).  The  squares  correspond  to  the  different  cohorts  in 
the  simulation.  They  give  the  fraction  of  the  total  population  in  a  cohort  of  mean  maturity 
stage  p  and  mean  cell  quota  q.  (w  =  1,  U  =  10,  [po?Pc]  =  [0.1, 0.4],  s;  =  5,  d  =  0.5) 
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Figure  2-16:  Population  dynamics  for  decreasing  dilution  rates.  For  d  <  0.2,  the  model 
converges  to  steady-state.  As  d  decreases,  the  oscillatory  transients  increase  in  amplitude. 
Eventually,  for  d  =  0.1  the  oscillations  persist,  (v  =  1,  U  =  10,  \po,Pc]  =  [0.1, 0.4],  Si  —  5). 
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Another  critical  parameter  is  the  inflowing  nutrient  concentration  S{.  As  s,-  increases 
the  steady-state  is  replaced  by  persistent  oscillations  of  increasing  amplitude.  This  is  shown 
in  Figure  2-18  where  the  amplitude  of  the  oscillation  in  total  cell  numbers  is  plotted  as  a 
function  of  s,-. 


Figure  2-18:  Bifurcation  diagrams  for  an  increasing  nutrient  inflow.  As  st-  increases  the  long¬ 
term  dynamics  of  the  system  changes  from  a  steady-state  to  a  limit  cycle.  This  is  shown 
by  plotting  the  amplitude  (A)  and  the  coefficient  of  variation  (B)  of  total  cell  numbers  for 
different  values  of  s;.  (v  =  1,  U  =  10,  d  =  0.1,  [po>Pc]  =  [0.1, 0.4].) 

The  above  oscillations  are  generation  cycles  resulting  from  the  interaction  between  re¬ 
source  dynamics  and  the  population  distribution  along  the  cell  cycle.  Figure  2-19  illustrates 
this  point  by  following  a  complete  cycle  for  total  cell  numbers,  cell  numbers  in  the  segment 
[PoiPch  ambient  nutrient  concentration,  and  the  mean  cell  quota  in  [po,Pc]-  As  the  total 
number  of  cells  increases  (r  ~  7,  upper  panel),  the  ambient  nutrient  concentration  decreases 
by  uptake  (middle  panel).  Then,  the  mean  cell  quota  of  cells  in  [po^Pc]  decreases  (bottom 
panel)  and  cells  move  slowly  through  this  nutrient-dependent  segment.  As  a  result,  cells 
accumulate  in  \po,pc\  (t  —  10,  upper  panel)  and  do  not  proceed  towards  division.  Then, 
total  cell  number  decreases,  uptake  goes  down,  and  ambient  nutrient  levels  go  up.  The 
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resulting  increase  of  mean  cell  quota  in  [po,Pc]  allows  cells  to  proceed  through  the  cycle. 
The  next  pulse  in  total  numbers  follows  (r  ~  15,  upper  panel).  This  completes  a  generation 
cycle.  Notice  that  the  number  of  cells  in  [po,Pc\  does  not  oscillate  in  synchrony  with  total 
cell  numbers.  These  changes  in  population  distribution  along  the  cell  cycle  indicate  that 
the  oscillations  differ  from  typical  consumer-resource  fluctuations. 


6.0 

co 

a5 

■Q  4.0 

E 

ZJ 

c 

=  2.0 
CD 

o 

0.0 

0 


c 

CD 


rs 

c 

-*—• 

a 

<n 

E 

< 


Time  (x) 


Time  (x) 


Figure  2-19:  Generation  cycles.  The  oscillatory  dynamics  of  the  model  involves  changes  in 
the  population  distribution  along  the  cell  cycle.  The  interplay  of  this  distribution  with  the 
consumer-resource  interaction  drives  the  oscillations  and  synchronizes  the  cell  population. 
See  text  for  a  complete  description  of  this  figure,  (v  =  1,  U  =  10,  s;  =  5,  d  =  0.1, 
[Po,Pc]  =  [0.1, 0.4]). 
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Figure  2-20:  Square-wave  perturbation.  After  turning  the  chemostat  flow  off  and  then  on, 
total  cell  numbers  display  oscillatory  transients  (see  text  for  an  explanation).  This  type 
of  perturbation  produces  short-term  fluctuations  of  large  amplitude.  In  (A),  d  =  0.1  and 
the  oscillations  persist.  In  (B),  d  =  0.2  and  the  oscillations  slowly  decay,  (v  =  1,  U  =  10, 
Si  =  5,  \po,Pc)  =  [0.1, 0.4]). 


Because  the  model  exhibits  long  transients,  the  time  for  the  build-up  of  generation 
cycles  with  large  amplitude  may  be  long  (see  Figure  2-16  for  d  =  0.1,  and  Figure  2-17). 
The  amplitude  of  the  short  term  fluctuations  depends,  however,  on  the  type  of  initial 
perturbation  (i.e.  on  initial  conditions).  One  type  of  perturbation  used  to  study  transients 
in  the  chemostat  consists  of  turning  the  flow  off  and  then  on  some  time  later  (Williams, 
1971).  In  the  model,  this  so-called  square-wave  perturbation  sets  initial  oscillations  of 
large  initial  amplitude.  This  happens  because  without  any  flow  through  the  chemostat,  the 
population  rapidly  increases  and  consumes  the  resource.  Then,  starvation  follows  and  forces 
the  synchronization  of  the  cells  in  the  nutrient-dependent  segment  of  the  cycle.  When  the 
flow  is  turned  on  again,  this  synchronization  persists  and  is  reinforced  by  the  mechanism 
described  in  the  previous  paragraph.  Figure  2-20  illustrates  this  point  for  two  different 
values  of  the  dilution  rate.  Notice  that  even  for  a  dilution  rate  (d  —  0.2)  that  eventually 
leads  to  an  equilibrium,  transient  oscillations  of  large  initial  amplitude  occur.  As  is  shown 
below,  these  oscillations  are  important  in  the  response  of  the  system  to  a  variable  nutrient 
supply. 
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2.3.4  Model  dynamics  under  a  variable  nutrient  supply 

As  before,  I  consider  a  periodic  nutrient  supply  5,(l  +  esin(Y1t))  of  period  T  and  amplitude 
e.  The  non-dimensional  equation  for  the  ambient  resource  becomes 

sj  Q  C 

—  =  d(si(  1  +  esinwr)  -  s)  -  U  — — ntot.  (2.29) 

dr  i  +  s 

where  the  non-dimensional  frequency  o>  =  (2ir)/(Tv0)  and  the  mean  nutrient  inflow  s,  = 
Si/Ku ■  The  other  equations  remain  unchanged. 

Simulations  for  different  frequencies  u  have  shown  that  the  model  is  capable  of  quasiperi- 
odic  dynamics.  Thus,  the  population  displays  aperiodic  behavior  with  variability  at  frequen¬ 
cies  other  than  that  of  the  forcing.  This  is  shown  in  Figure  2-21  for  u  =  1  and  d  =  0.1.  The 
initial  condition  for  this  simulation  is  chosen  as  the  end-point  of  the  square-wave  perturba¬ 
tion  experiment  in  Figure  2-20(A)  .  Thus,  the  population  is  at  r  =  0  partially  synchronized 
and  has  an  intrinsic  frequency  of  oscillation  given  by  the  generation  cycles.  For  a  periodic 
nutrient  inflow,  total  cell  numbers  exhibit  two  dominant  frequencies.  Only  one  of  these 
equals  the  frequency  of  the  forcing  (compare  Figures  2-21  (A)  and  (B)).  Another  example  is 
given  in  Figure  2-22  for  d  =  0.2  with  initial  conditions  set  by  the  end-point  of  the  simulation 
in  Figure  2-20(B)  .  Notice  that  for  these  parameter  values  and  a  constant  nutrient  supply, 
the  system  eventually  converges  to  an  equilibrium.  For  a  periodic  supply,  however,  the  long¬ 
term  behavior  of  total  cell  numbers  is  aperiodic  (Figure  2-22  (B)).  This  irregular  behavior 
is  also  apparent  in  the  short-term  dynamics  and  is  not  simply  tracking  the  nutrient  forcing 
(Figures  2-22(A)  and  (C)).  To  demonstrate  the  aperiodic  nature  of  the  dynamics,  one  of  the 
variables  is  plotted  against  itself  at  lagged  intervals  of  time.  The  resulting  trajectory  moves 
on  the  surface  of  a  torus  and  never  repeats  itself  (Figure  2-23).  The  corresponding  power 
spectrum  displays  variability  at  frequencies  other  than  that  of  the  forcing  (Figure  2-24  ). 
The  cross-correlation  between  population  patterns  and  nutrient  forcing  is  low  at  any  time 
lag  (Figure  2-25  ). 

2.4  Other  extensions 

Other  extensions  of  the  original  model  were  considered  to  investigate  the  robustness  of  its 
qualitative  behavior.  Uptake  limited  to  part  of  the  cycle  did  not  modify  the  main  results. 
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Time  (x) 


Figure  2-21:  Population  response  to  a  periodic  nutrient  forcing.  The  periodic  nutrient  sup¬ 
ply  is  shown  in  (B).  The  population  response  is  quasiperiodic.  The  long  period  modulation 
of  these  patterns  is  not  present  in  the  environmental  forcing,  (v  =  1,  U  =  10,  s;  =  5, 
d  =  0.1,  [p0,Pc]  =  [0.1, 0.4],  e  =  0.9,  w  =  1). 
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Figure  2-22:  Population  response  to  a  periodic  nutrient  forcing.  The  periodic  nutrient 
supply  is  shown  in  (C).  The  dynamics  of  total  cell  numbers  is  quasiperiodic  (A,B)  and 
displays  variability  at  frequencies  other  than  that  of  the  forcing.  The  irregular  population 
patterns  are  shown  for  both  the  transient  (A)  and  long  term  dynamics  (B).  (v  =  1,  U  =  10, 
=  5,  d  =  0.2,  [po, Pc]  =  [0.1, 0.4],  e  =  0.9,  w  =  1). 
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Figure  2-23:  Torus  attractor.  The  attractor  of  the  system  can  be  reconstructed  by  plotting 
one  of  the  variables,  the  ambient  nutrient  s,  against  itself  at  lagged  intervals  of  time  after 
transients  have  been  removed.  Three  dimensions  (s(r),  s(t  +  1.5)  and  s(t  +  3))  suffice  to 
reconstruct  the  attractor  since  trajectories  move  on  the  surface  of  a  torus,  (v  =  1,  U  =  10, 
si  =  5,  d  =  0.2,  [p0,  Pc]  =  [0.1, 0.4],  e  =  0.9,  u>  =  1). 


Frequency  (2 ni)' 

Figure  2-24:  Power  spectrum  of  quasiperiodic  population  dynamics.  Total  cell  numbers 
display  variability  at  frequencies  other  than  the  forcing  frequency  u>.  The  power  spectrum 
is  obtained  for  the  dynamics  of  total  cell  numbers  after  removal  of  transients  (v  =  1,  U  =  10, 
Si  =  5 ,d  =  0.2,  \po,pc]  =  [0.1, 0.4],  e  =  0.9,  w  =  1). 
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Similarly,  the  form  of  the  division  probability  distribution  <f>{p)  seemed  to  have  little  influ¬ 
ence  on  the  qualitative  behavior  of  the  model.  I  considered  normal  distributions  of  smaller 
variance  as  well  as  beta  distributions  for  4>(p).  In  the  latter  case,  all  cells  divide  before  p  =  1 
which  is  defined  as  the  maximum  maturity  stage  of  a  cell.  The  corresponding  division  rate 
B(p)  is  an  increasing  function  of  p  which  becomes  arbitrarily  large  as  p  approaches  1. 

2.5  Discussion 

This  work  demonstrates  that  the  life-history  structure  (the  stages  of  the  cell  cycle)  can 
introduce  a  wider  range  of  dynamics  than  that  of  unstructured  models  for  the  nonlinear 
interaction  between  a  phytoplankton  population  and  a  limiting  nutrient  resource.  In  the 
chemostat  models  studied  here,  oscillations  in  population  numbers  (transients  or  not)  are 
possible  under  a  constant  nutrient  supply.  These  generation  cycles  involve  changes  in  the 
population  distribution  along  the  cell  cycle;  changes  driven  by  the  interplay  of  this  distri¬ 
bution  with  environmental  resource  levels.  Under  a  periodic  nutrient  supply,  cell  numbers 
can  exhibit  aperiodic  behavior  with  variability  at  temporal  scales  different  from  that  of  the 
forcing.  This  complex  response  occurs  because  the  generation  cycles  introduce  a  temporal 
scale  intrinsic  to  the  population  capable  of  interacting  with  environmental  forcing  frequen¬ 
cies.  Oster  and  Takahashi  (1974)  found  generation  cycles  in  age-classified  models  for  insect 
parasite- host  systems  and  determined  with  a  linear  approach  the  importance  of  this  intrin¬ 
sic  cycle  to  dynamics  in  periodic  environments.  Oscillatory  behavior  caused  by  changes 
in  population  structure  has  also  been  found  in  an  age-classified  fishery  model  (Levin  and 
Goodyear,  1980),  in  a  model  for  the  interaction  of  Daphnia  and  phytoplankton  (De  Roos, 
1992),  and  more  recently,  in  a  model  for  an  age-structured  cell  population  and  its  limiting 
nutrient  resource  (Minkevich  and  Abramychev,  1994).  Here,  I  have  related  observations  on 
resource  control  of  cell  cycle  progression  to  the  occurrence  of  generation  cycles  in  phyto¬ 
plankton  dynamics. 

Chemostat  models  without  population  structure  that  consider  all  cells  as  equal,  such 
as  the  Monod  and  Droop  equations,  do  not  display  any  oscillatory  dynamics  (Lande  and 
Oyarzun,  1992;  Cunningam  and  Nisbet,  1980;  Smith  and  Waltman,  1994).  Furthermore,  the 
Droop  model  exhibits  a  simple  response  to  a  periodic  supply  of  nutrients:  it  oscillates  at  the 
forcing  frequency  (Pascual,  1994).  Thus,  there  is  no  transfer  of  variability  to  other  temporal 
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scales  and  the  population  basically  tracks  the  environmental  forcing.  It  is  worth  noting  that 
these  chemostat  models  were  derived  from  steady-state  observations,  and  were  originally 
developed  for  the  dynamics  of  total  biomass  although  they  have  been  used  extensively  in 
the  literature  for  total  cell  numbers.  Droop  (1979)  warned  against  applying  the  models  to 
cell  numbers  and  to  transient  dynamics. 

When  applied  to  cell  numbers,  however,  these  traditional  models  fail  to  account  for 
the  transient  oscillations  observed  in  chemostat  experiments  (Cunningham  and  Maas,  1978; 
Caperon,  1969;  Williams,  1971).  Hypotheses  to  explain  such  oscillations  resorted  to  the  idea 
of  a  time  delay  between  resource  levels  and  population  growth  (Caperon,  1969,  Cunningham 
and  Nisbet,  1980;  Williams,  1971  ).  This  study  suggests  a  link  between  observations  at  the 
cell  cycle  level  and  such  delays.  In  fact,  the  population  distribution  along  the  cell  cycle 
introduces  a  variable  time  delay  between  resource  levels  and  population  growth  by  division: 
the  number  of  cells  undergoing  division  becomes  a  function  of  the  cells’  past  resource  history. 
Learning  about  resource  control  of  cell  cycle  progression  appears  then  relevant  to  a  better 
understanding  of  the  postulated  delays  in  the  nutrient-phytoplankton  interaction.  The  cell 
cycle  may  also  provide  an  explanation  for  the  multi-peaked  responses  observed  in  chemostat 
experiments  with  a  pulsed  nutrient  supply  (Olson  and  Chisholm,  1983).  The  occurrence  of 

aperiodic  behavior  remains  to  be  determined. 

The  cell  cycle  models  in  this  work  present  dynamics  consistent  with  the  observation  of 
transient  oscillations  in  chemostat  experiments.  Two  aspects  of  their  dynamics,  however, 
do  not  match  the  observations.  First,  the  transient  oscillations  described  in  the  literature 
display  a  much  faster  damping  (Cunningham  and  Nisbet,  1980;  Williams,  1971).  Second, 
evidence  for  persistent  cycles  is  weak  and  relies  only  on  comments  about  failed  experiments 
(Droop,  1966;  Pickett,  1974).  It  is  possible  that  the  models  need  to  be  improved  and  that 
more  needs  to  be  known  about  resource  control  of  cell  cycle  progression  (see  below).  It  is 
also  possible  that  the  main  qualitative  results  obtained  with  the  models  would  be  found 
under  the  appropiate  experimental  conditions.  To  speculate  on  this  point,  it  is  interesting 
to  cite  two  comments  from  the  chemostat  literature: 

1.  ‘The  chemostat  was  generally  very  unstable.  The  degree  of  instability  is  indicated 
by  the  fact  that  Table  I  reports  measurements  for  146  out  of  538  days  of  operation. 
Much  of  the  unreported  time  was  spent  in  periods  of  oscillation  of  more  than  50%  in 
cell  density.  The  stability  might  have  been  increased  by  reducing  Sa  at  the  expense 
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of  reduced  cell  density’  (Pickett,  1975); 

2.  ‘The  steady-sate  was  normally  reached  within  10  days  of  altering  the  dilution  rate, 
provided  the  alteration  was  less  than  about  20%.  Greater  alterations  tended  to  set 
up  oscillations,  particularly  at  the  lower  dilution  rates.  At  rates  below  0.2  vol.  per 
day,  these  oscillations,  which  might  be  of  rather  large  amplitude  and  of  long  period, 
tended  to  persist  and  a  true  steady-state  was  difficult  to  achieve.’  (Droop,  1966). 

Notice  that  the  experimental  conditions  leading  to  these  ‘persistent’  cycles  are  exactly  the 
ones  predicted  by  the  models  in  this  study:  high  nutrient  supply  and  low  dilution  rates. 

One  important  avenue  in  improving  the  models  would  be  a  better  representation  of 
nutrient  uptake.  This  would  require  a  population  model  that  follows  not  only  the  position 
along  the  cell  cycle  but  the  size  of  an  individual  cell;  that  is,  an  understanding  of  how  the 
processes  of  development  (the  division  cycle)  and  growth  (in  size  or  biomass)  are  coupled 
in  phytoplankton  cells.  Cell  cycle  models  that  incorporate  this  coupling  do  exist  for  other 
eucaryotic  cells  (Tyson,  1985).  In  phytoplankton,  however,  the  nature  of  the  coupling 
remains  unclear  and  the  transition  point  may  not  involve  a  size  treshold  (Olson  et  al. , 
1986). 

Would  a  better  representation  of  uptake  modify  the  main  conclusions  of  this  work? 
The  uptake  feedback  from  the  population  to  the  resource  is  an  essential  component  of  the 
generation  cycles.  Recall  that  these  oscillations  occur  because  a  pulse  in  division  is  followed 
by  an  increase  in  uptake  and  the  consequent  decrease  in  resource  levels  slows  progression 
of  the  cells  through  the  cell  cycle.  Thus,  it  is  critical  that  an  increase  in  cell  numbers  lead 
to  a  decrease  in  the  resource.  In  the  models,  this  occurs  because  uptake  is  a  function  of 
cell  numbers.  Uptake  differences  among  the  cells  were  ignored.  If  uptake  were  proportional 
to  biomass  and  uptake  were  continuous  during  the  cell  cycle,  generation  cycles  would  not 
occur  because  a  pulse  in  division  would  not  change  total  biomass.  Any  of  the  following 
two  conditions,  however,  would  produce  an  increase  in  uptake  after  a  pulse  in  cell  division. 
The  first  one  consists  of  uptake  being  restricted  to  a  segment  of  the  cycle.  Brezinski  (1992) 
has  recently  shown  that  the  uptake  of  silica  in  a  diatom  species  is  restricted  to  Gl.  The 
second  condition  consists  of  uptake  being  proportional  to  cell  surface  area.  A  pulse  in 
division  would  clearly  increase  the  total  uptake  area  of  the  population.  In  cells  that  obtain 
nutrients  by  the  process  of  diffusion,  uptake  is  a  function  of  their  surface  area  (Reynolds, 


55 


1994  (Appendix  6.1.2)).  Under  either  of  these  conditions,  I  expect  the  qualitative  results  of 
the  models  to  hold. 

Other  open  areas  related  to  this  work,  relevant  to  population  patterns  in  the  field, 
include  the  modelling  of  multispecies  phytoplankton  assemblages,  and  the  interaction  of 
light  and  nutrients  in  driving  the  cell  cycle.  While  chemostat  models  regard  population 
dynamics,  oceanographic  models  are  concerned  with  multispecies  assemblages.  An  open 
question  is  how  to  aggregate  species  with  similar  life  histories  in  a  multispecies  model  that 
incorporates  information  at  the  cell  cycle  level.  Heath  and  Spencer  (1985)  developed  a 
cell  cycle  model  driven  by  the  photoperiod  for  a  diatom  species.  Heath  (1988)  then  used 
this  model  for  a  phytoplankton  assemblage  in  the  field  and  pointed  out  that  variability  in 
the  duration  of  cell  cycle  stages  within  a  species  may  be  comparable  to  variability  among 
species. 

In  phytoplankton,  internal  sources  of  temporal  regulation  in  division  patterns  include 
circadian  rythms  and  the  cell  cycle  under  environmental  control  of  nutrients  and  light 
(Prezelin,  1992).  The  importance  of  these  different  sources  varies  with  the  group  of  algae 
(Prezelin,  1992),  and  the  interaction  of  nutrients  and  light  in  driving  the  cycle  is  not  well 
understood  (Heath,  1988).  Light  is  a  well  known  driving  force  of  the  cell  cycle  in  the 
field,  one  that  generates  regular  patterns  in  division  (as  do  circadian  rythms).  Experiments 
have  shown  that  in  some  species,  nutrient  forcing  can  override  the  patterns  produced  by 
the  light-dark  cycle  so  that  fluctuations  in  population  growth  rate  become  phased  to  the 
nutrient  pulses  (Olson  and  Chisholm,  1983;  Putt  and  Prezelin,  1988).  For  such  a  species, 
my  work  suggests  that  cell  numbers  may  display  irregular  patterns  in  the  field,  and  that  it 
may  be  difficult  to  infer  from  those  patterns  the  scale,  or  even  the  nature,  of  the  underlying 
environmental  forcing. 

The  importance  of  within  population  heterogeneity  has  been  demonstrated  repeatedly 
for  a  variety  of  plants  and  animals  (Caswell,  1989).  Unicellular  algae  have  been  neglected 
because  measuring  the  heterogeneity  (the  stages  of  the  cell  cycle)  has  been  difficult  until 
the  recent  developments  of  flow  cytometry  (Chisholm  et  a/.,  1986).  This  work  indicates 
that  such  heterogeneity  may  introduce  an  important  temporal  scale  of  population  response 
to  environmental  variability. 
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2.6  Appendix 


2.6.1  The  EBT  for  the  basic  cell  cycle  model 

I  briefly  describe  here  the  application  of  the  escalator  boxcar  train  method  (EBT)  to  solve 
numerically  the  basic  cell  cycle  model.  For  a  general  description  of  the  method  see  De  Roos 
et  al.  (1992);  for  a  more  technical  description,  see  De  Roos  (1988)  or  De  Roos  et  al.  (1988). 

The  basic  idea  of  the  EBT  is  to  approximate  numerically  the  dynamics  of  a  structured 
population  by  a  set  of  ordinary  differential  equations  (ODEs),  which  can  be  subsequently 
integrated  with  one  of  the  many  well-known  methods  available  (here,  a  4th  order  Runge- 
Kutta  method  with  a  variable  stepsize).  For  the  approximation,  the  structured  population 
is  subdivided  into  groups  of  individuals  called  cohorts.  A  cohort  consists  of  individuals  born 
within  the  same  interval  of  time,  called  cohort  interval.  Every  cohort  is  characterized  by 
a  specific  set  of  statistics:  the  total  number  of  individuals  in  the  cohort,  and  the  mean  of 
the  state  variable(s)  used  to  classify  the  population  (for  example,  age,  size,  or  maturation 
stage).  Given  any  time  interval,  there  are  two  fundamentally  different  types  of  cohorts:  the 
cohort  in  creation  which  contains  the  individuals  born  during  that  interval,  and  the  internal 
cohorts  made  of  individuals  born  before  that  interval.  One  can  visualize  the  EBT  as  a  two 
step  process  that  repeats  in  time:  during  a  cohort  interval  the  ODEs  for  the  statistics  of 
the  cohorts  are  simulated;  at  the  end  of  a  cohort  interval,  a  cohort  in  creation  becomes  an 
internal  cohort.  An  interesting  property  of  the  EBT  is  the  dynamic  character  of  the  total 
number  of  cohorts:  cohorts  are  created  at  a  rate  determined  by  the  cohort  interval,  and 
cohorts  are  eliminated  if  they  become  empty  or  contain  a  negligible  number  of  individuals. 

Let  7 ij  and  pj  denote  respectively  the  number  and  the  mean  maturation  stage  of  cells 
in  cohort  j,  and  let  j  =  0  correspond  to  the  cohort  in  creation.  I  rewrite  the  system  of 
equations  2.11,  2.12,  2.13,  and  2.14  as  the  following  system  of  ODEs. 

For  the  internal  cohorts  (j  ^  0), 


drij 

dr 

dPj 

dr 

dPj 

dr 


—dtij  —  b(pj)rij 


s 

1  +  5 


if  Pj  e  [ Po,Pc ] 


(2.30) 

(2.31) 


v 
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otherwise 


For  the  cohort  in  creation  (j  =  0),  the  number  of  cells  decreases  due  to  ouflow  and 
increases  due  to  reproduction, 


=  -dno  +  2  ^  b(pj)rij. 


(2.32) 


Instead  of  following  the  dynamics  of  po,  an  equation  is  written  for  the  new  variable  txq  = 
noPo-  This  is  necessary  because  the  mean  quantity  po  is  undefined  when  the  cohort  in 
creation  is  empty  (i.e.  no  =  0)  (see  Appendix  in  De  Roos  et  al.  (1992)  for  a  detailed 
description  of  the  problem  posed  by  no  =  0  in  the  equations  for  mean  state  variables). 


Then, 


ro  , 

—  =  vno  —  UTT  0- 


(2.33) 


When  the  cohort  in  creation  becomes  an  internal  cohort,  7r0  is  divided  by  no- 
Finally,  the  dynamics  of  the  resource  is  given  by 


—  =  d(si  -  s) - -  y2  rij . 

dr  y  ’  s  +  l^  J 

3 


(2.34) 


2.6.2  The  EBT  for  the  cell  cycle  model  with  storage 

The  EBT  formulation  for  the  cell  cycle  model  with  storage  introduces  a  new  variable:  the 
mean  cell  quota  qj  in  cohort  j.  Then,  the  dynamics  of  the  internal  cohorts  is  given  by 


-drij  —  b(pj)tij, 


(2.35) 


0  if  Pj  G  [p0,Pc]  and  qj  <  1 

(1  -  — )  if  pi  e  \poiPc]  and  q  >  l 
v  otherwise, 


(2.36) 

(2.37) 


dqj_  =  s 
dr  1  +  5 


(2.38) 
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For  the  cohort  in  creation,  ODEs  are  written  for  n0,  7Tq  =  p0n0  and  tTq  =  qono, 


=  -dn0  +  2  Y  KPj)nj 


(2.39) 


=  VUq  -  dlTfy 


(2.40) 


=  -dr o  +  U +  Y  b(Pj)nj(lj- 


Finally,  the  equation  for  the  resource  dynamics  becomes 


dr  y  ’  1  +  s  4^ 


(2.41) 
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Chapter  3 


Periodic  response  to  periodic 
forcing  of  the  Droop  equations  for 
phytoplankton  growth' 


If  thou  (dear  reader)  art  wearied  with  this  tiresome  method  of  computation, 
have  pity  on  me,  who  had  to  go  through  it  seventy  times  at  least,  with 
an  immense  expenditure  of  time... 

—Johannes  Kepler,  1609.  Astronomia  Nova 


3.1  Introduction 

In  the  ocean,  the  microscopic  algae  or  phytoplankton  are  faced  with  a  highly  variable  supply 
of  their  essential  nutrients  (Harris,  1980;  Kilham  and  Hecky,  1988).  The  method  of  contin¬ 
uous  culture,  known  as  the  chemostat  (Tempest,  1970),  provides  an  experimental  system 
to  investigate  the  consequences  of  this  variability  for  population  dynamics.  Phytoplankton 
ecologists  view  the  chemostat  as  the  most  simple,  yet  controllable,  idealization  of  an  aquatic 
system  with  both  an  inflow  and  an  outflow  of  nutrients. 

Equations  modelling  phytoplankton  population  dynamics  in  a  chemostat  originally  re¬ 
lated  the  growth  rate  of  the  cells  to  the  nutrient  concentration  in  the  medium,  as  described 
by  Monod  (1942)  for  microorganisms.  Later,  Droop  (1968,  1973)  modified  this  relation  by 
proposing  that  nutrient  uptake  was  a  function  of  the  ambient  nutrient  concentration,  but 
growth  rate  varied  with  the  internal  nutrient  level  of  the  cells. 

Most  studies  of  these  models  have  focused  on  either  steady-state  growth  under  a  con- 

1This  chapter  was  published  in  J.  of  Math.  Biol.  (1994)  32:743-759. 
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stant  nutrient  flux  or  on  the  transient  response  to  a  single  perturbation  (Burmaster,  1978). 
Turpin  et  al  (1981)  studied  the  effect  of  nutrient  fluctuations  on  phytoplankton  growth,  but 
simplified  the  model  with  steady-state  assumptions.  In  this  work,  I  focus  on  periodic  nu¬ 
trient  fluctuations  and  investigate  their  consequences  for  population  dynamics  with  the  full 
nonlinear  model  proposed  by  Droop.  In  the  model,  nutrient  storage  by  the  cells  introduces 
time  delays  between  the  environmental  nutrient  pool  and  population  growth.  These  time 
delays  have  the  the  potential  to  interact  with  the  periodic  supply  of  nutrients  to  generate 
a  complex  population  response.  I  show  that  this  is  not  the  case:  the  population  oscillates 
on  the  same  frequency  as  the  nutrient  forcing.  The  existence  of  this  oscillatory  solution  is 
proven  by  closely  following  the  approach  of  Cushing  (1977),  Butler  and  Freedman  (1981), 
and  Bardi  (1981),  to  models  of  predator-prey  interactions  in  periodic  environments.  A  posi¬ 
tive  periodic  solution  is  shown  to  bifurcate  from  a  trivial  solution  that  loses  stability.  These 
two  cycles  are  shown  to  exchange  local  stability  at  the  bifurcation  point.  Numerical  results 
indicate  that  the  positive  cycle  attracts  all  positive  trajectories. 

This  work  establishes  a  basis  for  future  comparison  of  the  model  to  experimental  data. 
Some  related  but  unsolved  theoretical  questions  are  briefly  discussed. 

3.2  The  model 

In  a  chemostat,  nutrients  at  an  input  concentration  St,  are  supplied  by  a  through  flow  at 
rate  F  into  a  chamber  of  volume  V.  The  effluent  contains  both  medium  and  phytoplankton 
cells,  and  the  residence  time  of  the  cells  in  the  chamber  is  given  by  the  reciprocal  of  the 
dilution  rate  D  =  F/V.  The  chamber  is  assumed  to  be  well-mixed  although  in  practice 
organism  growth  on  the  chamber  walls  may  violate  this  assumption. 

Three  state  variables  describe  the  dynamics  within  the  chemostat  chamber:  the  phyto¬ 
plankton  biomass  concentration  X  (biomass  per  unit  volume),  the  concentration  of  limiting 
nutrient  S  (mass  per  unit  volume),  and  the  concentration  of  limiting  nutrient  in  the  inter¬ 
nal  pool  Q  (also  known  as  the  cell  quota,  in  mass  per  unit  biomass).  In  these  definitions, 
biomass  can  be  replaced  by  cell  density  only  if  the  average  mass  of  a  cell  remains  fairly 
constant  in  time  (Droop,  1979).  Phytoplankton  growth  rate,  proceeds  at  rate  p,  while 
nutrient  uptake  proceeds  at  rate  p.  The  phytoplankton  death  rate  is  assumed  negligible 
in  comparison  to  the  washout  rate.  The  following  equations  (Droop  1968,  1973),  model 
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phytoplankton  growth  in  a  chemostat 


dr 

dS_ 

dr 

dQ 

dr 


fiX  -  DX 
D(Si  -  S)-pX 
p-  pQ 


with 

p  ”  Mm(l  *  hq/Q'j  and  p  —  pm 

where  pm  denotes  the  maximum  uptake  rate,  Kq  the  minimum  cell  quota,  pm  the  maximum 
nutrient  uptake  rate,  and  Kp  the  half  saturation  constant. 

With  the  dimensionless  variables 


Kp  +  SJ 


x 0-R. 

Kn  ha 


s  =  — — ,  and.  t  =  pmr 

hP 


the  model  becomes 


x  = 

s  = 

q  = 


x(l - )  —  ux 

q 

u{si  —  s)  -  Ux 


-9+1 


(3-1) 


where  dot  denotes  the  time  derivative  with  respect  to  t.  Instead  of  six  parameters,  the 
dimensionless  equations  contain  the  three  parameters 


Both  the  dimensionless  dilution  rate  u  and  the  dimensionless  nutrient  inflow  Si  are  under 
experimental  control.  By  constrast,  U  characterizes  a  phytoplankton  species  with  respect  to 
a  particular  nutrient  by  comparing  the  temporal  scales  of  population  growth  and  nutrient 
uptake.  Experimental  work  with  various  phytoplankton  species  shows  this  parameter  to 
cover  a  broad  range  of  values  (U  ~  1  —  200)  (DiToro,  1980). 

The  phase  space  of  biological  relevance  (in  which  equations  3.1  describe  the  chemostat 
system)  is  given  by  (x  >  0,s  >  0,9  >  1).  It  is  positively  invariant:  Consider  an  initial 
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condition  in  this  phase  space.  From  the  first  equation  of  3.1,  x  -  z(0)  exp(/o(l  -  u  -  l/q)), 
and  therefore,  x  remains  positive  for  all  time.  Then,  from  3.1,  s  also  remains  positive  and 
q  remains  larger  than  one. 


3.3  Analysis  of  the  model 

3.3.1  A  constant  nutrient  supply 

To  motivate  the  analysis  for  periodic  s,-,  the  behavior  of  system  3.1  for  constant  s,-  is  briefly 
described  in  the  region  of  parameter  space  where  u,  and  U  are  positive,  and  where  the 
maximum  growth  rate  exceeds  the  dilution  rate,  that  is,  u  <  1.  Outside  this  region  the 
population  cannot  persist  in  the  chemostat  chamber. 

In  this  case,  two  equilibrium  solutions  exist: 

Pi  = 

P2  = 

Consider  u  as  a  bifurcation  parameter  and  let  uc  —  (s,{7)/(l  +  Si(U  +  1)).  When  u  <  uc , 
the  trivial  solution  is  unstable  and  P2,  locally  stable  and  positive.  At  the  critical  value 
u  =  uc,  Pi  and  P2  coincide  and  exchange  stability.  For  u  >  uc,  growth  and  nutrient 
uptake  proceed  too  slowly  to  balance  cell  losses,  the  trivial  equilibrium  Pi  becomes  locally 
stable  and  P2)  now  negative,  loses  stability.  (For  a  proof  of  this  result  see  the  Appendix  or 
Lange  and  Oyarzun  (1992).  Because  Lange  and  Oyarzun  (1992)  consider  a  different  non- 
dimensional  form  of  the  equations,  I  have  presented  a  local  stability  proof  in  the  Appendix). 

Imagine  for  a  moment  that  the  existence  of  a  positive  equilibrium  P2  were  unknown.  It 
could  be  inferred  from  the  bifurcation  of  the  trivial  equilibrium  Pi  as  u  passes  through  uc. 
In  the  following  section,  a  similar  idea  underlies  the  proof  that  a  nontrivial  solution  of  the 
same  frequency  as  the  forcing  does  exist  for  periodic  s,-.  The  solution  is  shown  to  bifurcate 
from  a  trivial  cycle  that  loses  stability. 

3.3.2  A  periodic  nutrient  supply 

As  before,  I  consider  the  region  of  parameter  space  given  by  U  and  st-  positive,  and  u 
between  zero  and  one.  The  forcing  function  S{(t)  belongs  to  the  space  B,  defined  as  the 


0,  S{,  1  +  U 

x,  s,  q]  = 


1  +  Si 
(1  -  u)(Si  -  s), 


u-u(u  +  iy  i  —  u\ 
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Banach  space  of  T-periodic  continuous  functions  under  the  norm  \Y\0  =  sup0<t<r  |y(t)| 
where  T  is  an  arbitrary,  but  fixed  period.  The  notation  B 3  is  used  for  the  product  space 
B  X  B  X  B  under  the  norm  \X,Y,  Z\0  =  \X\0  +  |F|0  +  \Z\a.  Also,  for  Y  in  B ,  the  average 
of  Y  is  defined  as  (Y)  =  (1  /T)  /0T  Y{t)dt. 

The  trivial  solution 

Theorems  1  and  2  state  some  needed  results  on  the  trivial  solution  of  3.1,  the  solution  with 
no  cells  in  the  system. 

Theorem  1.  When  x  =  0,  system  3.1  admits  a  T-periodic  solution.  This  solution, 
denoted  by  (0,s*,q*),  satisfies  s*(t )  >  0  and  q*(t )  >  1  for  all  t  >  0  . 

Proof:  When  x  =  0,  system  3.1  becomes 
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sponding  to  deviations  from  the  trivial  cycle.  System  3.1  becomes 


£i(l  —  u - ) 

V  x3  +  q*J 


-UX  2  —  Ux\ 


-x3  +  u 


.  +  *2 
1  5*  +  x2  +  1 
3:2 +  S*  _  3* 

X2  +  s*  +  1  1  +  5* 


System  3.2  can  be  written 


where  the  functions 


zi(l  -  u - -)  +  h(xux3) 

9 

Us* 

-ux2 - — — -®i  +  f2{xi,x2) 

s*  +  1 

- X 3  +  (g«  +  ^*2  +  Mx*) 


h{x  1,0:2)  = 

Mx  2)  = 


q*  x3  +  q* 

Ux\s*  rT  s*  +  x2 

- —  —  U  X 1 - — 

s*  +  1  5*  +  x2  +  1 

-Ux2  (  X2  +  s* 

(  o*  1  ^2  ^  t  o*  _L  1 


(s*  +  l) 


X2  +  5*  +  1  S*  +  1  ' 


contain  higher  order  terms  arbitrarily  close  to  the  trivial  solution  (0,0,0).  This  is  shown  by 
the  following  series  expansions,  valid  when  \x2(t)\  <  s*(t )  +  1  and  |x3(t)|  <  q*(t)  for  all  t, 

OO 

>  =  S(-ir+,(^ 

«**.**>  =  -0  £(-^'0^-0  t(-ir 


fl 

M**)  =  PD-  1)"+*(8.+2l)n4,- 


Thus,  if  the  linear  system 


xUl  —  u - ) 

u  q* 

Us* 

—  UX  2 - .1 

S*  +  1 
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u 

x3  =  -*3  + 

(s*  +  l)2 

is  (locally)  uniformly  asymptotically  stable,  the  same  is  true  of  (0,0,0)  for  the  nonlinear 
system  3.3  (Halanay,  1966). 

Let  g(t)  =  1  —  u  —  1  /q*(t)  and  write  g(t)  =  (g)  +  A g.  Then, 
x\ (t)  ~  xi(0)efo  =  xi(0)e^te^o 

But  e^o  belongs  to  B,  and  therefore,  when  (5)  is  negative,  xi  tends  exponentially  to 
zero  as  t  becomes  arbitrarily  large.  Then,  by  the  second  and  third  equations  in  3.5,  the 
same  is  true  for  ar2  and  S3,  and  hence  3.5  is  (uniformly)  asymptotically  stable.  Conversely, 
if  (9)  >  0,  then  3.5  has  solutions  starting  arbitrarily  close  to  (0,0,0)  for  which  Sj(t)  does 
not  approach  0.  It  follows  that  (0,0,0)  is  unstable  for  (5)  when  ( g )  >  0. 


Bifurcation  of  the  trivial  solution 

Now  consider  what  happens  when  (g)  >  0  and  the  trivial  solution  loses  stability.  The 
following  theorem  states  the  main  result  on  the  existence  and  local  stability  of  a  positive 
cycle  of  exactly  the  same  frequency  as  the  nutrient  forcing.  (Here,  positive  refers  to  the 
state  variables  remaining  positive  for  all  time). 

Theorem  3.  When  (1  —  u  -  1  /q*)  >  0  there  exists  a  positive  T-periodic  solution  of 
system  3.1.  This  solution  is  locally  asymptotically  stable  for  values  of  u  arbitrarily  close  to 
uc  satisfying  (1  —  uc  —  1/q*)  =  0. 

Notice  that  the  condition  (g)  >  0  could  be  stated  as  a  condition  on  u  if  the  values  uc 
satisfying  (g)  =  0,  were  known.  Denote  the  smallest  such  u  in  (0, 1)  by  ucm.  The  following 
facts  about  (1  -  1/q*}  establish  that  (g)  >  0  when  u  belongs  to  (0,ucm).  First,  (1  -  1/q*) 
is  a  continuous  function  of  u  in  (0, 1).  Second, 


lim  (1  - 

u—+  0 


U(sj) 

U  ( S{ )  +  ( Si )  +  1 


and  therefore 

0  <  lim  (1  — — )  <  1. 
u— o'  q*1 

Thus,  (1  —  1/q*)  >  u  (or  equivalently  (g)  >  0)  for  u  in  (0,wcm).  Figure  3-1  illustrates 
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this  point  for  a  sinusoidal  forcing  function.  The  curve  y(u )  =  (1  —  1/q*)  intersects  the 
diagonal  y(u )  =  u  at  u  =  uc,  and  for  u  <  u,c,  (1  —  1  /q*)  >  u.  Figure  3-1  also  illustrates  that 
this  curve  crosses  the  diagonal  at  a  single  point.  Equivalently  the  root  of  ( g )  =  0  is  unique 
(i.e.  uc  =  ucm  is  unique).  This  result,  obtained  numerically  in  an  extensive  exploration  of 
parameter  space,  supports  the  conjecture  that  ( g )  >  0  for  u  <  uc  and  ( g )  <  0  for  u  >  uc. 
Equivalently,  when  u  is  reduced  below  the  critical  value  uc  the  trivial  solution  loses  stability. 


Figure  3-1:  The  curves  y(u )  =  (1  -  1/q*)  (...)  are  shown  for  s,  =  1  + 
0.9sin(0.2t)  and  for  different  values  of  the  parameter  U  (from  top  to  bottom:  U  = 
200,50,40,30,20,10,9,8,7,6,5,4,3,2,1).  Each  curve  intersects  the  diagonal  y(u)  =  u  at  a 
unique  point. 


Theorems  2  and  3  show  an  exchange  of  local  stability  at  u  =  uc  similar  to  the  one 
described  for  a  constant  nutrient  forcing.  Theorem  3  states  the  local  stability  of  the  positive 
periodic  solution,  both  dynamically  local  in  the  sense  of  linearized  stability  and  local  near 
uc.  It  does  not  address  the  global  stability  of  the  solution.  However,  extensive  simulation  of 
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system  3.1  indicates  that  when  u  <  uc ,  all  positive  trajectories  converge  to  this  cycle  and, 
when  u  >  uc,  trajectories  converge  to  the  origin.  Figure  3-2  shows  some  numerical  results 
for  a  sinusoidal  nutrient  forcing.  In  this  case,  the  simulations  also  confirm  the  critical  value 
uc,  estimated  here  as  falling  between  0.65  and  0.7. 


(A) 


Figure  3-2:  Phase  portrait  of  the  deviations  from  the  trivial  cycle  (0,s*,  q*)  for  different 
values  of  the  parameter  u.  The  number  by  each  curve  corresponds  to  u.  The  forcing  function 
is  Si  =  1  +  0.9(sin(0.2t),  and  U  =  5.  As  u  approaches  the  critical  value  uc  &  0.7,  the  limit 
cycles  approach  the  origin.  (In  (A):  projection  x\  =  x  vs.  x2  —  s  -  s* .  In  (B):  projection 
x3  =  q  -  q*  vs.  x2  =  s  —  s*.) 

To  prove  theorem  3,  the  following  three  lemmas  are  needed.  The  first  one,  due  to 
Cushing  (1977)  and  extended  here  to  include  one  more  variable,  concerns  the  existence  of 
periodic  solutions  for  a  particular  3-dimensional  system  with  periodic  coefficients.  This 
lemma  is  used  to  write  system  3.1  as  an  operator  equation  to  which  results  from  bifurcation 
theory  apply.  Lemma  2,  a  local  bifurcation  result  due  to  Krasnoselskii  (1964),  is  then  used 
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to  show  that  arbitrarily  close  to  a  trivial  cycle,  a  periodic  solution  exists.  Finally,  this  local 
result,  valid  for  the  bifurcation  parameter  arbitrarily  close  to  a  critical  value,  is  extended  to 
a  larger  region  of  parameter  space  by  applying  a  global  bifurcation  result  due  to  Rabinowitz 
(1971)  and  stated  in  Lemma  3. 

A  FEW  LEMMAS:  Lemma  1  (Cushing,  1977). 

Let  aij  G  B,  i,j  =  1,2,3. 

•(A)  If  (an)  0  for  i  =  1,2,3,  then  the  linear  homogeneous  system 


2/1  =  «n2/i 

V2  =  a22y2  +  a2iyi 

(3.6) 

il3  =  a33y3  +  a32y2 

has  no  nontrivial  solution  in  B3 . 

functions  fi  in  B, 

In  this  case,  the  nonhomogeneous  system 

with  forcing 

Xi 

—  an^i  +  /1 

x 2 

=  a22x2  +  a2\Xi  +  f2 

(3.7) 

x3  =  a33x3  +  a32x2  +  h 

has  a  unique  solution  (x\,  x2,  x3)  in  B3 .  If  L  denotes  the  operator  from  B3  to  itself,  assigning 
to  each  set  of  forcing  functions  (f\,f2->fz)  a  solution  (xi,x2,x3)  of  3.7,  then,  L  is  linear 
and  compact.  Furthermore,  if  Li  denotes  the  operator,  from  B  to  B,  mapping  the  forcing 
function  f  to  the  solution  of  +  f,  then  the  operator  L  may  be  decomposed  as 

L(fuf2,h)  =  (xi,x2,x3) 

—  {L\f\,L2(a2\L\fi  +  f2),  L3(a32L2(a2iLifi  +  f2 )  +  /3)). 

•  (B)  //(an)  =  0  and  (a22)  ^  0  ^  (0,33),  then  3.6  has  exactly  one  independent  solution 
in  B3 . 

In  the  next  two  lemmas,  G(X,x)  denotes  a  one- parameter  family  of  continuous  compact 
operators  from  E  =  R  X  X  to  X,  where  A  is  a  real  Banach  space.  Furthermore,  G  = 
A L(x)  +  H(X,x)  with  L  linear  and  compact  and  H,  o(||x|[)  for  x  near  0  uniformly  on 
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bounded  A  intervals. 

Lemma  2(Krasnoselskii,  1964;  Rabinowitz,  1971).  If  Xc  is  a  characteristic  value  of  L  of 
odd  multiplicity,  then  (Ac,0)  is  a  bifurcation  point  of  the  equation  G(X,x)  =  x  with  respect 
to  the  curve  of  trivial  solutions. 

Let  (1  be  an  open  set  in  E  containing  (Ac,0),  and  C  be  the  set  of  nontrivial  solutions 
of  G(X,x)  =  x  in  E.  The  following  lemma  states  a  global  bifurcation  result  for  the  case  of 
non-globally  defined  operators  H . 

Lemma  3  (Rabinowitz,  1971;  Bardi,  1984).  Assume  that  Xc  has  multiplicity  1,  and  that 
H  is  defined  on  Cl,  H  is  independent  of  X,  and  H  is  continuously  (Frechet)  differentiable  in 
a  neighbourhood  of  ( Ac,0).  Then,  C  contains  two  connected  branches  of  solutions,  meeting 
at  (Ac,  0),  and  each  satisfying  one  of  the  following  alternatives.  Each  branch: 

(i)  is  unbounded  in  E,  or 

(ii)  meets  dCl,  the  boundary  of  Cl,  or 

(iii)  meets  (A,  0)  where  X  is  a  characteristic  value  of  L,  (X  ^  Xc). 

Proof  of  Theorem  3:  In  the  following  proof  of  theorem  3,  a  new  real  parameter  A  is 
introduced  in  system  3.1  and  chosen  as  a  bifurcation  parameter.  A  more  natural  choice 
would  appear  to  be  u.  However,  when  3.1  is  written  as  an  operator  equation,  its  linear  part 
L  depends  on  u  but  not  on  A,  (see  below).  Thus,  the  application  of  the  above  lemmas  is 
simplified  by  introducing  A  and  considering  the  following  system 


x  = 

s  — 

Q  = 


z(l  -  A - ) 

q 

u(si(t )  -  s)  —  U  x  — 
_q+1  +  U-±- 


(3.8) 


Note  that  3.8  and  3.1  coincide  for  A  =  u.  To  prove  theorem  3, 1  will  show  that  for  A  smaller 
than  a  critical  value,  including  the  desired  case  A  =  u  when  (1  —  —  u)  >  0,  system  3.8 
and  therefore  3.1,  has  a  nontrivial  T-periodic  solution. 

System  3.1  can  be  written 

X\  =  *!(1  -  A- ^-)  + /x(a:i,i3) 
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~ux 2  -  ,  ,  Xi  +  f2(x1,  X2) 
s*  +  1 

U  ,  ,  , 

“*3  +  (^  +  1)2a:2  +  HX  2) 


where,  as  before,  the  variables  =  x,  x2  =  s  —  s*  and  X3  =  q  —  q*  correspond  to  deviations 
from  the  trivial  solution,  and  the  functions  /,■  are  defined  in  3.4. 

From  theorem  1,  q*  >  1,  and  therefore,  (1  —  !/<?*}  ^  0.  But  then,  the  linear  system 


*■(*-?) 

Us* 

~ uy 2  “  ^TT2'1 
“ y 3  + 


(3.10) 


satisfies  the  assumptions  of  Lemma  1(A): 


=  — «  ^  0 
=  -1#0 


With  the  operator  A  of  this  Lemma,  system  3.9  can  be  equivalently  written  as  the  operator 


equation 


(£i,a:2, £3)  =  AL*(£i,£2,£3)  +  H(xi,x  2, £3) 


(3.11) 


where 


H(x !,£2,£3)  =  {Lifi,L2  y-  ~rp[Lifi  +  f2J  , 


(s*  -f  l)2 


i2  +/a)  +/s))- 


In  3.11,  L*  :  B3  — ►  X?3  is  linear  and  compact,  and  XT  :  A?3  — >■  A?3  is  continuous  and  compact 
since  Xj,  X2  and  A3  are  compact.  Also  H  is  of  order  higher  than  linear  near  (0,0,0).  Thus, 
Lemma  2  applies  to  equation  3.11,  and  bifurcation  occurs  at  the  nontrivial  solutions  of  the 


75 


linear  equation  (2/1,  S/2, 2/3 )  =  AL*(i/i,  y2,  jft),  or  equivalently,  given  the  definition  of  L*,  at 
the  nontrivial  solutions  of 


in 

=  yi(!-A-^) 

Us * 

(3.12) 

V2 

=  uy*  .  .  ,yi 

s*  +  1 

U 

h 

2/3  +  (S*  +  l)2  2/2 

From  Lemma  1(A,B),  3.12  has  a  nontrivial  solution  in  B 3  if  and  only  if 


A  =  Ac  =  (1-4>-  (3-13) 

Since  Ac  can  be  shown  to  have  multiplicity  1,  (see  Butler  and  Freedman,  1981),  bifurcation 
does  in  fact  occur  at  this  characteristic  value.  Then,  equation  3.11  admits  a  continuum  of 
nontrivial  solutions  in  R  X  B 3,  forming  two  branches  that  meet  at  (Ac,  0,0,0). 

Near  the  bifurcation  point,  the  set  of  nontrivial  solutions  is  investigated  with  the  fol¬ 
lowing  Lyapunov-Schmidt  small  parameter  expansions, 


A  —  Ac  +  Axc  +  A2(c)c 

Xi  =  xne  +  Xi2<?  +  Xi3(t,e)e2  (i  =  1,2,3)  (3.14) 


where  f  is  a  small  parameter  and  |A2(c)|  =0(|c|),  |a;,-3(i, e)|0  =0(|e|).  Substituting  these 
series  in  3.9  (with  the  functions  /,•  written  in  expanded  form)  and  equating  coefficients  of  e 
and  e2,  one  obtains 

in 
*21 

*31 

and 

*12  =  *12(1  -  Ac  -  ^)  -  xn(Ai  -  jrppO-  (3.18) 

Given  an  initial  condition  £n(0)  >  0,  equation  3.15  implies  that 


^ll(l  *) 

a 

(3.15) 

IJs* 

UX21  *  ,  !*11 

s*  +  1 

(3.16) 

U 

131  V+i  r121 

(3.17) 
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£ii(0  =  a;ii(0)eJo(1  Ac  l!q  ^  is  positive  for  all  positive  t.  Then,  from  equations  3.16  and 
3.17,  both  £21  and  2:3!  are  negative.  Also,  Ai  must  be  negative,  since  xi2  in  equation  3.18 
belongs  to  B  if  and  only  if  Aj  =  {x31/(q*)2)  (Halanay  (1966)  p.226,  or  Lemma  2  in  Cushing 
(1977).  This  lemma  states  that  for  a  in  B  and  (a)  =  0,  the  equation  x  =  ax  +  /  ,  /  e  B, 

0  a^ds)  =  0).  It  follows  that  arbitrarily  close  to 
the  bifurcation  point,  the  two  branches  of  nontrivial  solutions,  denoted  respectively  by  Ct 
and  Cx ,  satisfy 

Ct  -  {(^,xi,x2,x3)  e  R  X  B3  :  Xc  -  b0  <  X  <  Xc  for  some  bQ,  x\  >  0,  x2  <  0,  x3  <  0} 

Cx  =  {(A,  x\,  x2,  x3)  £  R  X  B3  :  Xc  <  X  <  Xc  +  b0  for  some  ba,  x1  <  0,  x2  >  0,  x3  >  0} 

Let  C£  correspond  to  Cx  when  this  set  is  defined  with  the  variables  x,s,q  instead  of 

xi,  x2,  x3.  That  is, 

C2  =  {(A,  x,s,q)  e  R  x  B3  :  Xc  -  b0  <  X  <  Xc  for  some  b0,  x  >  0,  s  <  s*,  q  <  q*}. 

To  determine  the  existence  of  a  T-periodic  solution  when  A  —  u  ,  (that  is,  when  system 
3.8  and  3.1  coincide),  the  extension  of  the  branch  C£  is  investigated  globally  in  fl,  the 
subset  of  R  X  B3  in  which  H(X,x,s,q)  is  defined.  More  specifically,  by  establishing  that 
T-periodic  solutions  of  3.8  exist  for  A  in  the  whole  interval  (0,  Ac),  the  desired  case,  A  =  u,  is 
captured  for  all  u  satisfying  (1  -  u  -  l/q*)  >  0.  This  idea  is  sketched  in  Fig.  3-3.  Note  that 
in  the  parameter  space  X/u ,  the  critical  value  Ac  corresponds  to  the  curve  A(u)  =  (1  —  1  /q*). 
When  (1  —  u  —  1/ q*)  >  0,  the  curve  A(u)  =  (1  —  l/q*)  is  above  the  diagonal  A(u)  =  u.  Thus, 
for  a  fix  u ,  the  interval  (0,  Ac),  contains  the  point  A  =  u. 

Because  H(X,x,s,q )  is  defined  for  x3  +  q*  ^  0  and  s*  +  x2  +  1  #  0  (see  equations 
3.4),  the  subset  is  chosen  as  ft  =  {R  x  B3  :  q  >  0  and  s  >  -1}.  Let  C+  denote  the 
extension  of  in  R  X  B3,  and  the  sets  A  and  T  denote  the  projections  of  C+,  onto  R  and 
B3,  respectively.  The  branch  C+  satisfies  one  of  the  three  alternatives  of  Lemma  3.  The 
third  alternative  is  impossible  since  by  Lemma  1,  there  does  not  exist  another  characteristic 
value  of  L.  The  following  facts  about  solutions  of  3.8  establish  that  the  second  alternative 
is  also  impossible. 

From  the  first  equation  of  3.8,  a:  =  x(0)e[io  where  h  =  1-X-l/q,  and  therefore,  x 
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Figure  3-3:  A  sketch  of  the  proof  in  parameter  space  X/u.  The  two  curves  A  =  u  and 
A  =  (1  —  1/q*)  intersect  for  u  =  uc.  Fix  up  for  any  value  of  u  smaller  than  uc.  The  local 
bifurcation  result  shows  the  existence  of  T-periodic  solutions  arbitrarily  close  to  points  on 
A  =  (1  -  1/q*)  (see  A  =  (Up,  Ac)).  This  local  result  is  then  extended  to  point  B  on  A  =  u 
for  which  systems  3.8  and  3.1  coincide. 
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remains  of  the  same  sign  for  all  t  >  0.  Assume  that  there  exists  a  solution  in  C*,  with  x  <  0. 
Since  this  branch  is  connected  and  contains  solutions  with  x  >  0  in  C£,  one  obtains  the 
contradiction  that  the  trivial  solution  (Ac,  0,  s*,  </*)  belongs  to  C+,.  Thus,  for  all  solutions 
in  C4.  x  is  positive,  and  from  3.8,  s  is  also  positive  and  q  is  larger  than  one.  Hence,  C+, 
does  not  meet  the  boundary  of  Cl  and  the  second  alternative  of  Lemma  3  does  not  apply  to 
this  branch. 

Finally,  the  first  alternative  of  Lemma  3  must  hold  and  either  A  or  T  are  unbounded 
(i.e.  C+,  contains  solutions  with  \x\0  or  |A|  arbitrarily  large).  It  is  shown  next  that  A  is 
unbounded  below  and  therefore  contains  the  whole  interval  (0,AC). 

The  set  A  is  bounded  above  by  Ac,  (assuming  otherwise  implies  that  there  exist  solutions 
in  C+,  with  x  <  0).  Assume  that  A  is  bounded  below.  From  3.8,  if  h  is  written  as  ( h )  +  A/i, 
then  x(t)  =  x(0)e^{efo  Ahd Since,  e^o  Ahd ^  belongs  to  B ,  then  x  belongs  to  B  if  and  only 
if  ( h )  =  0.  In  addition,  since  q  >  1,  there  exists  a  constant  M  such  that  |/t|0  <  M,  for  all 
solutions  in  C+ .  Thus,  there  exists  constants  N  and  P  such  that  fg  Ahd£  <  N  and  therefore 
\x\ o  <  P  for  all  solutions  in  C+ .  But  then,  T  is  bounded,  which  contradicts  Lemma  3.  It 
follows  that  A  must  be  unbounded  below  and  therefore  contains  the  whole  interval  (0,AC). 
In  particular,  there  exists  a  solution  in  C+  for  A  =  u.  Thus,  system  3.8,  (or  equivalently, 
system  3.1),  has  a  periodic  solution  in  B3.  This  completes  the  proof  of  the  existence  of  a 
positive  periodic  solution. 

To  prove  the  local  stability  of  this  solution  near  uc  the  proof  of  theorem  8  in  Cushing 
(1982)  is  closely  followed.  Local  stability  is  demonstrated  for  A  arbitrarily  close  to  Ac.  By 
applying  this  result  to  values  of  Ac  arbitrarily  close  to  uc  (see  Figure  3-3),  one  obtains  the 
desired  local  stability  for  u  near  uc. 

Let  N(p )  denote  and  open  ball  in  R  x  B3  of  radius  p  >  0  and  center  (Ac,  0,0,0), 
and  let  C+  denote  the  extension  of  Cf  in  R  x  B3  .  The  following  arguments  show  that 
the  solution  (*i,  x%,  x3)  of  system  3.9  is  locally  asymptotically  stable  for  (A,  X\,  X2,  x3)  £ 
C+nN(p)-(  Ac,  0,0,0). 

To  determine  the  stability  properties  of  the  branch  solution  C+  arbitrarily  close  to 
(Ac,  0,0,0),  system  3.9  is  linearized  at  C+.  This  linearization  yields, 

^  ^  ^  x3  +  q*^Vl  +  (x3  +  q*yV3 
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2/2 


(3.19) 


The  local  stability  of  the  branch  solution  C+  arbitrarily  close  to  (Ac,  0,0,0)  is  determined 
by  the  Floquet  exponents  of  system  3.19.  When  the  solution  ( \,xi,X2,x3 )  is  written  with 
the  small  parameter  expansions  in  3.14,  these  Floquet  exponents  are  also  functions  of  the 
parameter  e.  Notice  that  for  e  =  0,  system  3.19  becomes  3.12  (with  A  =  Ac).  Since  3.12  is 
a  block  triangular  system,  two  of  its  Floquet  exponents  are  those  of  the  reduced  system 


2/2  =  -uy2 

yz  =  -2/3  +  (,+V  (3-20) 

(Cushing,  1982).  Since  3.20  is  locally  asymptotically  stable  at  (0,0)  ,  these  two  Floquet 
exponents  must  have  negative  real  parts.  Thus,  for  e  sufficiently  small,  two  Floquet  ex¬ 
ponents  of  3.19  must  also  have  negative  real  parts.  The  remaining  exponent  of  3.12  is 
((1  —  Ac  -  —))  =  0.  Thus,  one  needs  to  determine  the  location  in  the  complex  plane  of  the 
remaining  exponent  of  3.19  when  e  is  small. 

But  e  is  a  Floquet  exponent  of  3.19  if  and  only  if  the  system 


z\ 

h 

zz 


(1  - A  - e)zi + 

-«*2  -  et2  -  U -  11  -  U  A +:*1  .*1 

(5*  +  X2  +  l)  S*-fX2  +  l 


-z3  -  ez3  + 


(s*  +  x2  +  l)2 


(3.21) 


has  a  nontrivial  T-periodic  solution  for  Z{  €  B ,  (i  =  1,2,3)  (Cushing,  1982).  The  sign  of  the 
real  part  of  e  is  obtained  by  expanding  e  =  e*e  +  e2(e)e  and  Zj  =  zn(t)  +  z,-2(t)e  +  Zj3(t,e)e 
where  e  is  a  small  parameter  and  where  |ei(e)|  =0(|e|)  and  |z,-3(t,  e)|0  =0(je|).  Substituting 
these  series  in  3.21  and  equating  coefficients  for  the  lowest  order  terms,  one  obtains  equations 
for  zn(t)  equivalent  to  equations  3.15,  3.16  and  3.17  for  xn(t),  (i  =  1,2,3).  Thus,  zn(t)  is 
positive  for  all  t  and  both  z2i(t)  and  z3i(t)  are  negative  for  all  t.  Equating  coefficients  for 
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the  first  order  e  terms  yields 


il2  -  (1  -  Ac  -  — )*12  +  (-Al  +  -  eX)zn  +  ^2^3!- 


Then,  212(0  belongs  to  B  if  and  only  if 


/  .  ,  *31  *31*11  ,  _  n 

+  7777?  _  ei  _  .  /..Vi)  “  °> 


(?*) 


zn{q*¥ 


(Lemma  2  in  Cushing  (1977)),  or  equivalently, 


ei  =  ( 


*31*11 

*n(<r)2 


) 


where  the  value  of  Ai  has  been  used.  Thus,  e\  <  0  and  e(c)  <  0  for  solutions  in  C+  fl  N(p). 
It  follows  that  solutions  in  C+  fl  N(p)  -  (Ac,  0,0,0)  are  locally  (uniformly  asymptotically 
)  stable.  For  Ac  arbitrarily  close  to  uc,  C+  fl  N(p)  contains  solutions  of  system  3.2  with 
((1  —  u  —  ^-))  >  0  (see  Figure  3-3).  This  completes  the  proof  on  the  local  stability  of  the 
T-periodic  solution  of  system  3.1. 


3.4  Discussion 

In  spite  of  its  nonlinearity,  the  Droop  model  predicts  a  simple  response  of  a  phytoplankton 
population  to  the  periodic  supply  of  nutrients:  an  oscillation  of  exactly  the  same  frequency 
as  the  forcings.  The  above  results  prove  the  existence  and  local,  but  not  global,  stability  of 
this  periodic  solution.  However,  extensive  simulation  of  the  model  indicates  that  positive 
trajectories  converge  to  this  cycle  in  the  same  region  of  parameter  space  where  the  trivial 
solution  is  known  to  be  unstable.  Also,  the  population  response  appears  to  mimic  the 
temporal  pattern  of  the  nutrient  inflow  (Fig.  3-4). 

This  simple  response  of  the  model  to  a  periodic  nutrient  supply  may  relate  to  its  also 
simple  behavior  under  constant  forcing.  In  fact,  the  Droop  equations  under  a  constant 
supply  of  the  resource  present  no  oscillations  in  their  approach  to  equilibrium  (Lange  and 
Oyarzun,  1992).  Thus,  the  system  lacks  any  internal  frequency  capable  of  interacting  with 
external  fluctuations  to  generate  complex  dynamics. 

These  results  provide  a  basis  to  evaluate  the  model  against  experimental  data  in  studies 
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(A) 


Figure  3-4:  Numerical  solutions  for  two  different  forcing  functions.  Curves  for  the  nutrient 
inflow  Si  and  the  phytoplankton  biomass  x  are  shown.  (st-  corresponds  in  (A)  to  the  sine 
function  1  +  0.95  sin(0.4f),  and  in  (B),  to  the  pulse  function  £)*l!f(l/(10(£  —  24 k)2  +  0.01)). 
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of  phytoplankton  population  dynamics  and  nutrient  variability.  I  am  aware  of  a  single  such 
study  by  Olson  and  Chisholm  (1983).  Their  data  for  Hymenomonas  carterae  grown  under 
daily  pulses  of  ammonium,  reveals  the  expected  presence  of  a  24  hour  periodicity.  However, 
a  simple  cycle  with  the  same  temporal  pattern  as  the  forcing  may  not  be  the  whole  story. 
After  filtering  and  averaging  the  data,  Olson  and  Chisholm  (1983)  observe  more  than  one 
single  peak  within  24  hours.  More  experiments  are  needed  to  provide  longer  data  sets  and 
cover  a  larger  range  of  forcing  frequencies. 

Finally,  two  open  questions  related  to  this  work  are  briefly  mentioned.  The  first  one 
concerns  extending  the  results  to  general  functions  of  uptake  and  growth.  This  generaliza¬ 
tion  would  include  extensions  of  the  Droop  model  that  consider  nutrient  uptake  a  function 
of  both  the  ambient  and  cellular  nutrient  levels  (Rhee,  1973).  The  second  one  regards  a  dif¬ 
ferent  approach  to  phytoplankton  dynamics  that  views  the  population  as  distributed  along 
the  cell  cycle.  Since  nutrient  levels  are  known  to  affect  progression  of  a  phytoplankton  cell 
through  its  cycle  (Vaulot  et  al  1987),  this  distribution  may  play  an  important  role  in  the 
population  response  to  nutrient  fluctuations.  This  approach  would  consider  the  dynamics 
of  both  population  biomass  and  cell  numbers  in  fluctuating  environments. 

Models  of  the  interplay  between  environmental  variability  and  phytoplankton  dynamics 
in  continuous  culture  benefit  from  their  powerful  experimental  setting.  Their  significance 
extends,  however,  to  the  broader  context  of  oceanographic  models  that  incorporate  the 
nutrient-phytoplankton  interaction.  It  is  therefore  essential  that  chemostat  models  capture 
the  transfer  of  variability  from  the  environment  to  the  population. 
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Appendix: 

Local  stability  of  the  trivial  steady-state  Pi:  The  Jacobian  matrix  of  model  3.1 


u  0  0 


J  =  ~uvk 


-U  0 


has  eigenvalues  given  by  its  diagonal  elements.  These  eigenvalues  are  all  negative  if  and 


only  if 


>  SiU/(si(U+  1)+  1). 


Local  stability  of  the  steady-state  P2'  Let  Nt0t(t)  =  xq  +  s  denote  the  total  con¬ 
centration  of  nutrients  in  the  culture  at  time  t.  Then,  from  3.1, 


=  -uNtot  +  usi 


and  Ntot  =  Ce  ut  +  S;  for  some  constant  C.  As  t  — ►  00,  Ntot  —■ >  Si  and,  provided  that  x  ^  0, 
system  (1)  becomes  equivalent  to  the  two  dimensional  system 


x(l - - — )  —  ux 

S{  —  s 

u{si  -  s)  -  Ux — — 


(3.22) 


At  steady-state,  the  Jacobian  matrix  of  3.22 


u  —  1  —(1  -  u )2 

_u_  |  (U-uU-u)2,  _ u _ \ 

u— 1  '  u—  1  1  U(l—u)—u' 


has  a  positive  determinant  if  and  only  if  S{  >  yr~-  But  then,  since  the  trace  of  J  is 
negative,  the  steady-state  is  locally  stable. 
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Chapter  4 

Diffusion-induced  chaos  in  a 
spatial  predator-prey  system- 


You  can’t  know  how  happy  I  am  that  we  met, 

I’m  strangely  attracted  to  you. 

— Cole  Porter.  It’s  All  Right  with  Me. 


4.1  Introduction 

A  variety  of  ecological  models  exhibit  chaotic  dynamics  because  of  nonlinearities  in  popu¬ 
lation  growth  and  interspecific  interactions  (e.g.  Gilpin,  1979;  Hastings  and  Powell,  1991; 
Kot  et  al.,  1992;  Schaffer,  1988).  These  models  have  for  the  most  part  ignored  space.  Ex¬ 
plicit  consideration  of  space,  however,  can  fundamentally  alter  the  dynamics  of  nonlinear 
interactions  (Turing,  1952;  Levin  and  Segel,  1976;  Segel  and  Jackson,  1972). 

The  few  ecological  studies  of  chaos  in  spatial  systems  consider  models  in  discrete  time 
and  space  (Sole  and  Vails,  1992;  Hassell  et  al.,  1991)  or  in  discrete  time  and  continuous  space 
(Kot,  1989).  In  all  these  models,  the  diffusive  dispersal  of  organisms  drives  a  predator-prey 
or  a  host-parasitoid  system  into  chaotic  dynamics. 

The  results  of  discrete  models  cannot  be  applied  directly  to  nonlinear  interactions  and 
dispersal  in  continuous  time  and  space.  It  is  well-known  that  discrete  models  exhibit  chaos 
more  readily  than  their  continuous  counterparts.  For  instance,  chaotic  dynamics  is  possible 
for  discrete  time  models  of  even  a  single  species,  but  require  at  least  three  variables  in 
continuous  time. 

1This  chapter  was  published  in  Proc.  R.  Soc.  Lond.  B  (1993),  251:  1-7. 
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In  this  paper  I  investigate  the  behavior  of  a  continuous  predator-prey  system  on  a  spa¬ 
tial  gradient  that  affects  the  intrinsic  growth  rate  of  the  prey.  This  model  differs  from 
most  reaction-diffusion  models  because  space  is  heterogeneous  rather  than  homogeneous. 
It  differs  from  coupled  map  lattices  and  other  discrete  models  because  it  is  continuous  in 
both  time  and  space.  Unlike  the  discrete-time  models  underlying  coupled  map  lattices  , 
the  predator-prey  models  used  here  cannot  exhibit  chaos  in  the  absence  of  spatial  diffu¬ 
sion.  Thus  any  chaotic  behavior  must  result  from  the  interaction  of  the  (non-chaotic)  local 
dynamics  with  the  spatial  gradient. 

Simulations  of  the  model  have  indicated  that  diffusion  can  drive  predator  and  prey  num¬ 
bers  into  complex  patterns  of  variability  in  time.  The  main  goal  of  my  work  is  to  determine 
if  these  patterns  are  chaotic.  I  demonstrate  with  a  variety  of  approaches,  including  bi¬ 
furcation  diagrams,  correlation  dimension  estimates  and  Poincare  sections  of  reconstructed 
attractors,  temporal  power  spectra,  and  dominant  Lyapunov  exponents,  that  predator  and 
prey  numbers  at  a  fixed  spatial  location  exhibit  temporal  chaos  and  quasiperiodicity.  At 
some  spatial  scales  these  results  may  apply  to  planktonic  organisms  transported  by  turbu¬ 
lent  diffusion. 


4.2  The  model 

To  pose  the  problem  in  its  simplest  form  consider  a  single  spatial  dimension  along  which 
both  species  diffuse  at  the  same  constant  rate  D.  At  any  point  X  and  time  T ,  the  dynamics 
of  the  prey  (P(X,T))  and  predator  ( H(X,T ))  populations  are  given  by  a  reaction- diffusion 
model  with  logistic  growth  of  the  prey  and  a  type  II  functional  response  of  the  predator: 


dP 

dT 

8H 

dT 


R,P{  1  - 


CiP 
C2  +  P 


H - MH +D 


ACiP  „02P 

c7TPH  +  DdT2 

d2H 


dx 2' 


(4.1) 


The  parameters  Rx,  K,  M  and  1/A  denote  the  intrinsic  growth  rate  and  carrying  capacity 
of  the  prey,  the  death  rate  of  the  predator  and  the  yield  coefficient  of  prey  to  predator, 
respectively.  The  constants  C\  and  C2  parameterize  the  saturating  functional  response. 

To  describe  an  environment  surrounded  by  dispersal  barriers,  I  assume  zero  flux  at  the 
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boundaries.  Hence,  at  X  =  0  and  X  =  L, 


dP  _  dH 
dX~dX 


for  all  T. 


(4.2) 


In  the  absence  of  diffusion,  (4.1)  is  a  standard  predator-prey  system,  which  exhibits 
stable  equilibria  or  limit  cycles  (May,  1973). 

A  simple  form  of  environmental  heterogeneity  can  be  introduced  by  allowing  the  pa¬ 
rameters  in  (4.1)  to  vary  with  X .  In  this  chapter  I  consider  the  case  where  the  prey  rate  of 
increase  Rx  is  a  linear  function  of  X . 

This  chapter  is  concerned  with  the  effects  of  diffusion  and  heterogeneity  on  a  system 
which,  in  the  absence  of  those  factors,  exhibits  limit  cycle  dynamics.  There  is  a  large 
literature  on  the  related  problem  of  the  diffusive  instability  of  fixed  points  (Turing,  1952; 
Murray,  1989;  see  Levin  and  Segel  1976  for  a  predator-prey  example).  It  is  worth  noting 
here  that  the  conditions  for  such  instabilities  are  not  satisfied  by  system  4.1  (see  Discussion). 

The  model  can  be  simplified  by  introducing  the  dimensionless  variables  p  =  P/K  and 
h  =  AH/K.  Space  is  scaled  by  the  total  length  of  the  gradient  L,  and  time  is  scaled 
by  a  characteristic  value  of  the  prey  growth  rate  R.  Thus,  x  —  XjL  and  t  =  RT  where 
R  =  Rx{X 0)  for  some  Xq  in  (0,  L ).  System  4.1  becomes 


dp 

dt 

dh 

dt 


rxp(l-p)  - 


ap  L  ,  jd^p 

1  +  bp  +  dx* 


-^P-h-mh  + 

1  +  bp 


P— 

d8x 2 


(4.3) 


where  the  new  parameters  are 


Rx  ,  r  C\K  K  M  J  D 

rT  =  -=-  =  e  +  tx,  a  =  ■  — ,  b  -  — ,  m  =  -=r,  a  =  — -=. 

x  R  J  ’  C2R  C2  R  L2R 


At  the  boundaries,  now  given  by  x  =  0  and  x  =  1, 


dp  _  dh 
dx  dx 


for  all  t. 


(4.4) 


(4.5) 
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4.3  Numerical  methods 


Although  the  dynamics  of  (4.1)  in  the  absence  of  diffusion  are  well-understood,  there  is 
little  analytic  theory  for  the  system  in  space.  The  equations  were  integrated  numerically 
with  an  implicit  scheme  using  100  spatial  grid  sites.  This  scheme  combines  a  fully  implicit 
method  for  the  diffusion  terms  (Roache,  1972),  with  a  fourth-order  Runge-Kutta  method 
for  the  predator-prey  interaction  terms. 

For  nonlinear  equations,  the  discretization  introduced  by  numerical  methods  may  gener¬ 
ate  spurious  results.  To  test  this  possibility,  the  resolution  of  the  simulations  was  increased 
in  space  and  time,  and  the  system  was  integrated  with  a  different  numerical  scheme,  a  finite 
difference  method  (forward  in  time  and  centered  in  space).  In  all  cases  the  same  qualitative 
results  were  obtained. 

In  the  absence  of  diffusion,  the  simulations  match  the  known  behavior  of  the  system, 
i.e.  stable  equilibria  or  limit  cycles.  The  accuracy  of  these  periodic  solutions  was  used 
as  a  criterion  for  periodicity.  If  successive  maxima  coincided  to  the  fourth  decimal  place, 
solutions  were  considered  periodic. 

4.4  Results 

The  results  presented  here  are  based  on  numerical  analysis  for  a  set  of  parameters,  ( a  =  5, 
b  =  5,  to  =  0.6,  e  =  2  and  /  =  -1.4),  chosen  to  obtain  limit  cycles  at  each  fixed  location 
along  the  gradient  in  the  absence  of  diffusion.  One  diffusion  rate,  ( d  =  10~4),  was  initially 
studied.  Figure  4-1  illustrates  the  irregular  temporal  and  spatial  behavior  of  prey  numbers 
after  transients  have  died  out.  Predator  and  prey  densities  at  any  fixed  location  in  space 
are  aperiodic  in  time  (Figure  4-2)  and  show  sensitivity  to  initial  conditions  (Figure  4-3). 
The  following  results  focus  on  characterizing  these  irregular  motions,  determining  if  they 
are  quasiperiodic  or  chaotic,  and  documenting  the  bifurcations  produced  by  changes  in  the 
diffusion  rate. 

Transitions  to  chaotic  behavior  are  known  to  occur  along  different  routes  as  a  parameter 
is  varied  (Schaffer,  1988;  Schuster,  1984).  Identifying  one  of  these  known  routes  provides  a 
diagnostic  for  chaos.  The  diffusion  rate  d  was  varied  in  a  range  that  covered  a  rich  range 
of  dynamics  including  not  only  the  irregular  behavior  described  above  but  also  periodic 
solutions.  System  (4.3)  exhibits  periodicity  at  both  low  and  large  values  of  d,  ( d  ~  10~8 
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Figure  4-1:  Complex  spatiotemporal  patterns  in  prey  density.  Transients  have  been  re¬ 
moved,  (d  =  10-4). 


Time 

Figure  4-2:  Irregular  temporal  behavior  of  predator  and  prey  densities  at  location  x  =  0.85 
(d=  10-4). 


0.8 


Figure  4-3:  Sensitivity  to  initial  conditions.  Temporal  solutions  diverge  for  small  initial 
differences.  Trajectories  of  prey  numbers  at  x  =  0.85  are  shown  for  two  initial  conditions 
differing  by  0.001  in  both  p  and  h  at  every  x,  (d  =  10~4). 

and  d  ~  3  X  10-3,  respectively).  To  obtain  a  bifurcation  diagram  successive  local  maxima 
at  a  fixed  location  were  plotted  as  a  function  of  the  corresponding  diffusion  rate  for  1(T4  < 
d  <  3  x  10-3  (Figure  4-4).  A  period-1  trajectory  produces  a  single  point.  More  generally, 
periodic  trajectories  produce  finite  number  of  points.  Successive  maxima  of  quasipenodic 
and  chaotic  trajectories  spread  over  a  range  of  values.  Whereas  the  former  densely  cover 
this  range,  the  latter  present  a  complex  structure.  Note  the  different  qualitative  regions  m 
the  diagram.  For  large  values  of  the  diffusion  rate,  dynamics  are  periodic  (Figure  4-4,  a). 
For  smaller  values  of  d,  periodicity  is  lost  and  the  maxima  visit  a  whole  segment  (Figure 
4-4,  a,b).  At  two  points  in  the  diagram,  the  maxima  suddenly  scatter  on  a  larger  segment 
(see  arrows  in  Figure  4-4,  c,d).  Where  windows  appear,  trajectories  become  periodic  again 
(see  Figure  4-4,  c,  for  d  «  5  x  10~4).  The  windows  in  the  diagram  were  investigated  at 
a  higher  resolution  in  d.  Results  (not  shown  here)  reveal  clearly  periodic  behavior  within 
these  windows  over  a  range  of  d  values.  Figure  4-4  is  thus  reminiscent  of  the  quasiperiodic 
route  to  chaos.  In  this  route,  a  bifurcation  occurs  by  adding  a  second  frequency  to  a  periodic 
motion.  The  attractor  of  the  system  becomes  a  two  dimensional  manifold,  the  surface  of  a 
torus.  When  the  ratio  of  these  two  fequencies  is  rational  the  trajectory  on  this  surface  closes 
after  a  finite  number  of  cycles.  This  periodic  motion  is  called  a  frequency-locked  state.  For 
an  irrational  ratio  the  motion  is  called  quasiperiodic,  the  trajectory  never  closes  and  covers 
the  whole  torus.  After  quasiperiodicity,  a  transition  to  chaos  and  the  break  up  of  the  torus 
into  a  strange  attractor  becomes  possible  (Schuster,  1984). 
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Figure  4-4:  Bifurcation  diagrams.  Successive  maxima  of  prey  density  at  x  —  0.85  are  plotted 
for  increasing  values  of  the  diffusion  rate  after  transients  have  died  out.  (In  (a),  values  of  d 
differ  by  10~5,  in  (b),  by  IQ-6.) 
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DIFFUSION  RATE 
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Figure  4-5:  Bifurcation  diagrams  (cont.).  Successive  maxima  of  prey  density  at  x  —  0.85' 
are  plotted  for  increasing  values  of  the  diffusion  rate  after  transients  have  died  out.  (Values 
of  d  differ  by  10-6) 
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Visualizing  the  system’s  attractor  would  permit  this  hypothetical  scenario  to  be  inves¬ 
tigated.  The  infinite  dimensionality  of  system  (4.3)  precludes  a  simple  plot  of  trajectories 
in  phase  space.  However,  if  the  attractor  itself  is  low  dimensional,  one  may  reconstruct  it 
from  knowledge  of  a  single  variable  (Takens  1981;  for  ecological  discussion  see  Kot  et  al. , 
1988).  Suppose  that  the  attractor  lies  in  an  n-dimensional  space,  but  that  one  follows  only 
the  dynamics  of  a  single  variable  z(t).  Then,  for  almost  every  time  lag  r,  the  attractor  of 
the  ^-dimensional  time  series 

Z(t)  =  [z(t),z(t  +  r),z(t  +  2 r),...,z(t  +  (E  -  l)r)] 

is  qualitatively  the  same  as  the  unknown  attractor  of  the  n-dimensional  system  (Takens, 
1981;  Kot  et  al.  1988).  The  ’embedding  dimension’  K,  which  needs  to  be  sufficiently  high 
but  not  larger  than  2 n  +  1,  corresponds  to  the  notion  of  degrees  of  freedom,  in  the  sense 
of  providing  a  sufficient  number  of  variables  to  specify  a  point  on  the  attractor  (Farmer, 
1982). 

Theoretical  and  numerical  results  have  indicated  that  the  attractors  of  many  infinite 
dimensional  dynamical  systems  are  of  finite  dimension  (Farmer,  1982).  To  explore  this 
possibility  for  system  (4.3),  the  fractal  dimension  of  the  attractor  was  estimated  by  com¬ 
puting  its  correlation  dimension  Dc  (Grassberger  and  Procaccia,  1983).  (A  description  of 
the  algorithm  can  be  found  in  Bingham  and  Kot,  1989).  Figure  4-6  shows  log-log  plots 
of  the  correlation  sum  vs.  length  scale.  These  curves  exhibit  linear  regions  with  slopes 
that  provide  an  estimate  of  the  correlation  dimension  Dc.  This  quantity  was  computed  for 
increasing  values  of  E  until  convergence  occurred.  For  d  =  10-3,  the  estimated  correlation 
dimension  converged  to  a  value  of  2.0  (characteristic  of  motion  on  a  torus)  for  E  >  3.  Thus, 
three  dimensions  appear  to  be  the  minimum  number  needed  to  reconstruct  the  attractor. 
For  d  —  10-4,  which  produced  the  complex  dynamics  of  Figures  4- 1-4-3,  the  attractor’s 
dimension  converged  to  Dc  =  3.2  for  E  >  7.  The  fractional  dimension  is  characteristic  of 
strange  attractors.  For  d  =  10-4,  then,  E  =  3  is  certainly  too  low,  but  may  still  provide 
information  on  qualitative  changes  in  dynamics.  The  attractor  was  reconstructed  in  three 
dimensions  from  the  time  series  of  prey  densities  at  the  fixed  location  x  =  0.85  for  decreas¬ 
ing  values  of  d.  After  the  loss  of  periodicity,  the  reconstructed  trajectories  appear  to  move 
on  a  torus  (Figure  4-7,  a).  In  the  windows  of  the  bifurcation  diagram,  frequency-locking 
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Figure  4-7:  Reconstruction  of  the  attractor.  The  attractor  is  reconstructed  in  three  dimen¬ 
sions  from  trajectory  of  prey  numbers  at  x  =  0.85.  (  In  (a),  aperiodic  motion  on  a  torus, 
( d  =  10“3).  In  (b),  periodic  motion  on  a  torus,  ( d  =  5.06  X  10~4).  In  (c),  aperiodic  motion 
on  the  projection  of  a  strange  attractor,  ( d  =  10~4)). 

occurs  (Figure  4-7,  b).  As  d  decreases,  the  reconstruction,  now  only  a  projection  of  the 
attractor,  becomes  more  complex.  (Figure  4-7,  c). 

Cutting  the  reconstructed  trajectory  with  a  plane  and  following  those  points  in  the  orbit 
that  intersect  the  plane  yields  a  Poincare  section  which  transforms  a  continuous  flow  into  a 
discrete  map  of  lower  dimension  (Kot  et  al.,  1988).  The  Poincare  section  of  a  torus  is  a  one 
dimensional  curve  topologically  equivalent  to  a  circle.  When  frequency-locking  occurs,  the 
Poincare  section  reduces  to  a  finite  set  of  points.  Figures  4-8(a,b)  show  these  phenomena 
for  the  attractors  of  Figures  4-7(a,b).  As  d  decreases,  the  torus  breaks  up  into  a  more 
complex  attractor,  (Figure  4-8,  c).  This  occurs  at  the  values  of  d  where  the  maxima  scatter 
into  larger  regions  in  the  bifurcation  diagram.  Given  the  low  embedding  dimension,  this 
qualitative  change  could  correspond  to  chaos  or  motion  on  a  a  higher  dimensional  torus. 
The  former  is  supported  by  the  non-integer  value  of  the  correlation  dimension. 
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Figure  4-8:  Poincare  sections.  Sections  are  shown  for  the  reconstructed  attractors  of  Figure 
4-6.  (  In  (a),  quasiperiodicity;  in  (b),  frequency-locking;  in  (c),  chaos). 

Further  evidence  for  chaos  is  provided  by  the  power  spectra  of  the  system.  The  power 
spectra  of  the  quasiperiodic  behavior  at  a  fixed  location  exhibits  sharp  peaks  corresponding 
to  the  fundamental  frequencies,  their  harmonics,  and  sums  and  differences  of  harmonics 
(Figure  4-9,  a).  Notice  the  increase  in  background  noise  levels  when  d  decreases  (Figure 
4-9,  b).  Such  continuous,  broadband  spectra  are  characteristic  of  chaos.  Suggestive  as  these 
results  may  be,  many  authors  consider  sensitivity  to  initial  conditions  as  the  crucial  defining 
feature  of  chaos.  The  sensitivity  apparent  in  Figure  4-3  was  quantified  by  calculating  the 
dominant  Lyapunov  exponent  Ai,  which  measures  the  long-term  average  rate  of  divergence 
of  nearby  trajectories.  Positive  values  of  Ai  indicate  sensitive  dependence  on  initial  con¬ 
ditions,  whereas  Xi  =  0  for  quasiperiodic  or  periodic  motion  and  <  0  for  stable  fixed 
points.  The  dominant  Lyapunov  exponent  Ai  was  computed  from  the  model  by  a  standard 
technique  for  a  system  of  differential  equations  (see  Wolf  et  al.  (1985)  for  a  description  of  the 
method  and  an  application  to  a  system  of  ordinary  differential  equations).  For  d  =  10-3, 
the  estimated  exponent  is  Aa  =  0.  For  d  =  NT4,  the  estimated  exponent  converges  to 
2.9  X  10~3  bits  per  unit  time  (or  Ai  =  4.9  X  10-2  bits  per  period  of  oscillation  at  x  =  0.85). 
This  confirms  the  apparent  quasiperiodicity  at  d  =  10-3  and  chaos  at  d  =  10-4. 
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Figure  4-9:  Power  spectral  density  of  prey  time  series  at  x  =  0.85.  In  (a),  quasiperiodicity 
(d  =  10-3).  In  (b),  chaotic  behavior  ( d  =  10~4). 

4.5  Discussion 

Diffusion  on  a  spatial  gradient  may  drive  an  otherwise  periodic  predator-prey  system  into 
quasiperiodic  or  chaotic  behavior.  In  the  absence  of  explicit  space  such  a  two-species  system 
is  only  capable  of  equilibria  and  limit  cycles.  Thus,  diffusion  and  spatial  heterogeneity 
introduce  qualitatively  new  types  of  behavior  for  this  predator-prey  interaction.  Previous 
ecological  examples  of  chaos  in  spatial  systems  consider  discrete  maps  for  one  or  two  species, 
a  class  of  models  already  capable  of  complex  dynamics  with  no  spatial  dimension  (Kot,  1989; 
Sole  and  Vails,  1992). 

The  spatial  distributions  of  predator  and  prey  in  this  model  are  not  simple  reflections  of 
the  underlying  environmental  gradient.  Because  of  nonlinearities,  environmental  variability 
is  transfered  to  other  spatial  scales  (see  Figure  4-1).  The  resulting  pattern  reflects  the  non- 
separable  effects  of  the  environment  and  the  biology.  An  attempt  to  classify  the  pattern  as 
autonomous  vs.  induced,  or  physical  vs.  biological  would  clearly  fail. 

These  results  differ  in  several  ways  from  previous  studies  of  pattern  formation  with 
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reaction-diffusion  equations  in  ecology  (Murray,  1989).  Those  studies  have  shown  that 
coupling  diffusion  to  nonlinear  ecological  interactions  can  generate  spatial  pattern  (Okubo, 
1980).  When  inhibiting  factors  (e.g.,  predators)  disperse  faster  than  do  activating  factors 
(e.g.,  prey),  diffusion  can  drive  the  system  into  a  new  asymptotic  state  that  is  non-uniform 
in  space  (Levin  and  Segel,  1976;  Segel  and  Jackson,  1972).  In  all  ecological  examples  of 
such  diffusive  instabilities,  the  resulting  spatial  pattern  has  been  either  constant  or  periodic 
in  time  (Levin  and  Segel,  1985;  Kishimoto  et  al.,  1983).  By  contrast,  system  (4.1)  produces 
chaotic  dynamics  and  spatial  patterns  that  are  continuously  changing  and  exhibit  long-term 
unpredictability. 

In  addition,  two-species  predator-prey  systems  exhibit  diffusive  instabilities  only  under 
special  conditions  (autocatalytic  prey  growth  rates  (Levin  and  Segel,  1976),  Allee  effects 
(Mimura  and  Murray,  1978),  or  density- dependent  predator  death  rate  (Okubo,  1980)). 
None  of  these  conditions  is  needed  for  the  generation  of  spatial  pattern  in  this  model. 
System  (4.3)  describes  the  predator-prey  interaction  with  a  standard  set  of  terms  commonly 
used  in  ecology  (May,  1973). 

Many  studies  of  predator-prey  or  host-parasitoid  systems  in  heterogeneous  environments 
have  concluded  that  dispersal  is  a  stabilizing  influence,  one  that  moderates  temporal  fluc¬ 
tuations  (Hastings,  1982;  McMurtie,  1978;  Taylor,  1990).  The  spatial  coupling  of  local 
fluctuations  may  give  rise  to  stable  equilibria  or  reduce  the  amplitude  of  the  oscillations 
at  the  local  and  regional  levels  (Comins  and  Blatt,  1974;  Crowley,  1981;  McLaughin  and 
Roughgarden,  1991).  In  this  study,  in  contrast,  local  oscillations  give  rise  to  complex  tem¬ 
poral  dynamics. 

The  scaled  diffusion  coefficient  d  is  a  critical  parameter  for  this  qualitative  change  in 
dynamics.  Chaos  and  quasiperiodicity  occur  for  d  on  the  order  of  10~4  to  10-3.  These 
orders  of  magnitude  are  plausible  for  turbulent  diffusion  in  aquatic  environments.  For 
example,  with  a  characteristic  growth  rate  R  =  10~5s_1  (or  one  division  per  day,  typical  of 
phytoplankton  growth),  the  diffusion  rate  D  must  satisfy  D/L 2  ^  10~8  or  10-9,  to  produce 
a  scaled  diffusion  rate  d  on  the  order  of  1(T3  to  10-4.  In  the  horizontal  dimension,  such 
values  of  D  occur  at  spatial  scales  of  10  to  100  km  (see  Figure  2.4  in  Okubo,  1980).  In  the 
vertical  dimension,  they  become  possible  at  scales  of  10  to  50  m  for  the  lowest  diffusivity 
values  estimated  in  the  mixed  layer  (Denman  and  Gargett,  1983).  Thus,  the  results  of  this 
paper  may  apply  to  planktonic  organisms  transported  by  turbulence. 
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Figure  4-10:  Dynamics  of  mean  prey  density.  The  mean  is  computed  over  the  whole  gradient 
after  transients  have  died  out,  ( d  —  10~4). 

These  results  suggest  two  directions  for  further  research.  The  first  is  the  interplay  be¬ 
tween  the  spatial  sampling  scale  and  the  detection  of  complex  dynamics  in  natural  environ¬ 
ments.  In  the  model,  predator  and  prey  numbers  exhibit  sensitivity  to  initial  conditions.  It 
remains  to  be  determined  if  this  fundamental  property  of  chaos  is  inherited  when  sampling 
the  system  at  larger  spatial  scales.  Mean  densities  over  the  entire  gradient,  for  instance, 
show  the  same  irregular  fluctuations  observed  at  fixed  locations  (Figure  4-10,  cf.  Figure 
4-2).  Furthermore,  an  understanding  of  the  relation  between  local  and  regional  dynamics 
in  this  system  should  be  relevant  to  discussions  on  chaos  and  persistence,  and  on  dispersal 
and  persistence  in  heterogeneous  environments. 

A  second  direction  concerns  the  critical  ingredients  for  generating  complex  dynamics 
in  the  model.  Preliminary  results  indicate  that  the  steepness  of  the  gradient  is  important. 
Increasing  the  slope  from  zero  produces  a  series  of  bifurcations  (not  shown  here)  suggestive 
of  a  quasiperiodic  route  to  chaos.  Other  mechanisms  of  spatial  coupling  besides  diffusion 
could  induce  temporal  chaos.  Evans  (1977)  briefly  commented  on  the  possibility  of  chaos  for 
a  planktonic  food  web  model  in  which  organisms  move  by  a  vertical  current  shear  combined 
with  migration  and  turbulent  mixing.  Models  with  other  forms  of  spatial  variability  and 
other  terms  to  describe  ecological  interactions  may  also  exhibit  chaos.  In  fact,  a  variety  of 
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nonlinear  dynamical  systems  are  capable  of  complex  dynamics  when  extended  into  space 
(Crutchfield  and  Kaneko,  1987;  Nicolis  and  Gaspard,  1990;  Vastano  et  al.  1990).  The 
ecological  conditions  for  which  space  would  induce  chaos,  however,  remain  to  be  explored. 
This  work  has  identified  as  particular  conditions  the  weak  coupling  by  diffusion  of  predator- 
prey  cycles  in  heterogeneous  space.  In  heterogeneous  environments,  spatially  induced  chaos 
may  play  an  important  role  in  the  generation  of  complex  spatio-temporal  patterns. 
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Chapter  5 


Spatial  pattern  in  a  predator-prey 
system  with  complex  dynamics 


5.1  Introduction 

Environmental  heterogeneity  has  long  been  recognized  as  a  principal  determinant  of  eco¬ 
logical  pattern  (Andrewartha  and  Birch,  1954;  Gleason,  1926).  This  is  particularly  true 
for  marine  systems  in  which  plankton  patterns  are  often  explained  by  the  variability  of 
the  physical  environment  (Mackas  et  al. ,  198  ;  Steele  and  Henderson,  1994).  In  terrestrial 
ecology,  environmental  fluctuations  have  also  been  at  the  center  of  long  standing  questions, 
such  as  the  origin  of  patchiness  and  the  importance  of  density-dependence  (Andrewartha 
and  Birch,  1954;  Deutschman  et  al .,  1993;  Hassell  et  al .,  1989). 

A  common  approach  to  identify  the  environmental  cause  of  ecological  pattern  consists 
of  matching  dominant  scales  of  variability.  A  dominant  scale  (frequency-1  or  wavelength) 
in  the  fluctuations  of  a  biological  variable  would  result  from  the  forcing  by  the  physical 
environment  at  a  similar  scale.  Cross-correlation  and  cross-spectral  analysis  are  examples 
of  extensively  used  methods  to  infer  causality  from  variability  at  similar  scales.  In  a  review 
of  the  literature  on  physical  processes  and  planktonic  ecosystems,  Denman  and  Powell 
(1984)  give  numerous  examples  of  successful  results  with  the  scale  matching  approach;  they 
point  out,  however,  that  as  often  as  not,  an  ecological  response  could  not  be  identified 
for  a  particular  physical  forcing.  One  possible  explanation  is  nonlinearity  (Denman,  1994; 
Denman  and  Powell,  1984).  Only  in  linear  systems  the  scales  of  the  response  typically  match 
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the  scales  of  the  forcings.  In  nonlinear  systems,  by  contrast,  variability  can  be  transfered 
to  widely  different  scales. 

There  is  ample  evidence  for  nonlinearity  in  population  growth,  ecological  interactions 
and  the  response  of  ecosystems  to  perturbations  (Dwyer  and  Perez,  1983;  Dwyer  et  al., 
1978;  Ellner  and  Turchin,  1995;  Turchin  and  Taylor,  1992).  Temporal  models  of  nonlinear 
ecological  interactions  have  demonstrated  that  population  fluctuations  may  occur  at  scales 
other  than  those  of  external  forcings.  For  instance,  it  is  well  known  that  temporal  predator- 
prey  models  under  periodic  forcing  are  capable  of  chaotic  dynamics,  that  is,  of  a  response  in 
which  all  temporal  scales  are  present  (Doveri  et  al.,  1993,  ;  Kot  et  al.,  1992;  Rinaldi  et  al., 
1993;  Schaffer,  1988).  In  fact,  chaos  displays  a  continuous  power  spectrum  with  variance  at 
all  frequencies.  This  complex  behavior  results  from  the  temporal  interplay  of  two  different 
oscillators:  the  predator-prey  cycles  and  the  seasonal  forcing  (Kot  et  al.  ,  1992). 

This  chapter  investigates  the  interplay  of  predator-prey  cycles  and  heterogeneous  space. 
It  develops  the  study,  initiated  in  Chapter  4,  of  predator-prey  interactions  and  diffusion 
along  environmental  gradients.  In  that  chapter,  I  showed  that  weak  diffusion  along  a  spatial 
gradient  may  drive  an  otherwise  periodic  system  into  complex  temporal  dynamics  that 
include  chaotic  behavior.  Here,  I  focus  on  the  spatial  properties  of  the  gradient  and  their 
consequences  for  the  spatio-temporal  dynamics  of  the  system.  In  particular,  I  ask:  how  do 
the  spatial  patterns  of  the  populations  compare  to  the  underlying  environmental  gradient? 
Throughout  this  chapter,  complex  dynamics  refers  to  chaos  and  quasiperiodicity. 

The  spatial  predator-prey  model  in  Chapter  4  is  formulated  as  a  continuous  model 
that  treats  space  as  an  explicit  and  continuous  variable.  More  specifically,  the  model  is 
a  reaction-diffusion  equation,  in  which  time  and  population  densities  are  also  continuous. 
Reaction-diffusion  equations  have  been  used  extensively  in  ecological  studies  to  investigate 
the  problem  of  pattern  formation  (see  Murray  (1989)  or  Okubo  (1980)  for  reviews).  The 
emphasis  has  been  on  understanding  how  ecological  interactions  would  lead  to  pattern  in  a 
spatially  homogeneous  environment  (Levin  and  Segel,  1976;  Segel  and  Jackson,  1972).  Few 
theoretical  studies  with  reaction-diffusion  equations  have  incorporated  heterogeneous  space 
by  considering  spatially-varying  parameters  or  diffusion  coefficients  (Benson  et  al.,  1993; 
Pacala  and  Roughgarden,  1982;  Pascual,  1993;  Malchow,  1993;  McLaughin  and  Roughgar- 
den,  1991,  1992).  Of  these,  only  the  model  of  Chapter  4  exhibits  spatio-temporal  chaos  or 
quasiperiodicity  (Pascual,  1993). 
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There  are  other  ecological  models  that  do  exhibit  complex  dynamics  in  space-time  (Has¬ 
sell  et  al. ,  1991;  Kot,  1989;  Sole  and  Vails,  1992).  Examples  comprise  models  of  predator- 
prey  and  host-parasitoid  interactions  that  treat  time  as  discrete.  These  discrete  models  differ 
from  the  model  studied  here  in  two  fundamental  ways.  First,  they  are  already  capable  of 
complex  dynamics  in  the  temporal  dimension:  space  and  dispersal  are  not  a  requirement  for 
chaos.  Second,  space  is  homogeneous.  The  predator-prey  model  used  here  cannot  exhibit 
complex  dynamics  in  the  absence  of  diffusion.  Chaos  results  from  the  interaction  of  the 
(non-chaotic)  local  dynamics  with  the  spatial  gradient  (Pascual,  1993).  Because  space  is 
heterogeneous,  the  model  provides  an  ideal  tool  to  investigate  how  the  population  patterns 
manifest  the  underlying  gradient  for  different  dynamic  regimes. 

This  Chapter  is  organized  as  follows.  First,  I  characterize  the  spatial  patterns  of  the 
populations  for  different  diffusion  coefficients.  Based  on  these  characterizations,  I  compare 
the  population  patterns  to  the  underlying  gradient  for  different  dynamic  regimes  (periodic¬ 
ity,  quasiperiodicity  and  chaos).  Then,  I  investigate  how  different  properties  of  the  gradient 
affect  the  spatio-temporal  dynamics  of  the  system.  This  identifies  properties  of  the  gradient 
that  are  essential  for  complex  dynamics.  I  finally  discuss  how  my  results  compare  to  other 
studies  of  predator-prey  interactions  and  diffusion  in  heterogeneous  space,  and  comment  on 
potential  implications  of  these  results  for  plankton  patterns. 

5.2  Spatial  pattern:  consequences  of  complex  dynamics 

5.2.1  Comparison  of  population  patterns  to  environmental  gradient 

In  Chapter  4,  I  showed  that  the  predator-prey  system  4.3  exhibits  a  sequence  of  qualitative 
changes  in  temporal  dynamics  as  the  diffusion  rate  d  decreases.  This  series  of  bifurcations 
corresponds  to  the  so-called  quasiperiodic  route  to  chaos:  as  d  decreases,  periodic  oscil¬ 
lations  give  place  to  quasiperiodic  fluctuations;  for  even  lower  values  of  d,  predator  and 
prey  dynamics  become  chaotic  (Pascual,  1993).  At  any  instant,  the  system  displays  spatial 
patterns  that  repeat  in  time  for  periodic  dynamics  but  continuously  change  for  chaos  and 
quasiperiodicity.  Figure  5-1  illustrates  typical  spatial  patterns  of  predator  and  prey  for 
two  different  diffusion  coefficients.  These  examples  are  suggestive  of  changes  in  scale,  as 
well  as  changes  in  the  degree  of  similarity  to  the  underlying  gradient,  with  different  dy¬ 
namic  regimes.  To  investigate  this  possibility  I  compared  the  population  patterns  to  the 
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underlying  gradient  with  the  methods  described  below. 


Figure  5-1:  Spatial  patterns  of  predator  and  prey  are  shown  for  chaotic  dynamics  (top 
panels,  d  =  10~4)  and  quasiperiodicity  (bottom  panels,  d  =  10-3)  at  three  different  times. 


Methods 

The  model  is  simulated  with  the  same  numerical  scheme  than  that  of  Chapter  4.  A  grid  of 
100  points  is  used  in  space.  Model  simulations  at  three  different  diffusion  coefficients  were 
selected  for  the  analysis  of  spatial  patterns: 

•  d  =  10-4:  The  system  exhibits  chaos.  Spatio-temporal  solutions  are  denoted  by  hc 
and  pc  for  predator  and  prey,  respectively. 

•  d  =  10-3:  The  system  exhibits  quasiperiodic  dynamics.  Solutions  are  denoted  by  hq 
and  pq  for  predator  and  prey,  respectively. 
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•  d  =  10  2:  The  system  exhibits  periodic  dynamics.  Solutions  are  denoted  by  hp  and 
pp  for  predator  and  prey,  respectively. 

All  other  parameters  in  these  simulations,  including  the  linear  gradient  in  the  prey’s 
growth  rate,  are  the  ones  of  Chapter  4.  The  above  simulations,  after  transients  are  removed, 
span  500  time  units  and  are  sampled  at  intervals  of  1. 

Two  additional  simulations  were  used  to  compare  cases  with  the  same  diffusion  coeffi¬ 
cient  ( d  =  10-4)  but  different  dynamic  regimes: 

•  A  gentler  gradient:  r(x)  =  —0.1a:  +  0.7  .  The  system  dynamics  are  periodic.  Solutions 
are  denoted  by  hgg  and  pgg  for  predator  and  prey,  respectively. 

•  A  higher  predator  mortality:  d  =  0.68.  The  system  converges  to  an  equilibrium 
solution  denoted  by  he  and  pe  for  predator  and  prey,  respectively. 


The  populations  patterns  were  compared  to  the  underlying  gradient  by  computing  a 
spatial  cross-correlation  coefficient  at  each  time  unit.  The  cross- correlation  coefficient  R 
between  two  variables  yi  and  x,-  (i  =  1  is  given  by 


where  the  overbar  denotes  the  mean.  Because  patterns  change  in  time,  the  correlation 
coefficients  also  vary.  I  therefore  plotted  a  histogram  of  these  coefficients  describing  the 
temporal  distribution  of  correlation  coefficients. 

Comparisons  were  also  made  by  obtaining  a  characteristic  scale  of  the  spatial  patterns  at 
each  time  unit,  and  then  plotting  a  histogram  of  characteristic  scales.  Here,  a  characteristic 
scale  is  defined  as  the  spatial  distance  one  has  to  travel  to  see  a  significant  change  in  the 
autocorrelation  function.  The  autocorrelation  function  measures  the  correlation  between 
values  at  different  distances  apart.  The  distance  at  which  the  autocorrelation  function 
crosses  0  gives  a  measure  of  characteristic  scale,  known  as  correlation  length  and  denoted 
here  by  Lc.  Let  yj  (j  =  1,2,. ..A)  be  a  data  set  with  values  sampled  at  a  distance  d.  Then, 
the  correlation  coefficient  for  values  separated  by  k  observations  (i.e.  a  lag  kd  apart)  is 
computed  as 


109 


(5.2) 


and  the  variance  Co  given  by 

<*  =  j  -  j)!  (5.3) 

3- 1 

(Chatfield,  1975).  The  correlation  length  Lc  =  dk  where  k  satisfies  r*,  =  0.  Note  that 
the  calculation  of  for  values  of  k  greater  than  about  N/4  is  questionable,  particularly 
when  the  goal  is  to  define  a  characteristic  scale.  This  is  because  the  pattern  contains  a 
trend  when  the  correlation  length  is  larger  than  a  fourth  of  the  total  extent  of  the  data 
set.  To  see  this,  consider  a  periodic  pattern  with  period  dN  (the  extent  of  the  whole  data 
set).  Then,  the  correlation  length  Lc  -  dN/ 4.  Patterns  with  longer  periods  do  not  repeat 
within  the  observed  space  frame  and  therefore  appear  as  a  trend.  For  a  non-stationary 
data  set,  a  trend  dominates  all  features  in  the  autocorrelation  function  and  the  values  of 
rjt  do  not  come  down  to  zero  except  for  large  values  of  the  lag  (Chatfield,  1975).  In  fact, 
the  sample  autocorrelation  function  is  considered  meaningful  only  for  stationary  data  sets. 
Here,  I  take  two  different  approaches  to  this  problem.  First,  for  the  purpose  of  comparing 
population  patterns  to  the  environmental  gradient,  which  is  clearly  non-stationary,  I  ignore 
the  stationarity  condition.  I  interpret  values  of  Lc  below  dN / 4  as  meaningful  estimates  of 
characteristic  scale  and  values  above  dN/4  as  indicators  of  a  trend.  The  longer  the  trend, 
the  more  the  pattern  resembles  the  linear  gradient.  In  fact,  for  samples  of  a  line,  r-jt  does 
not  cross  0  for  any  lag  in  (0,Nd).  Second,  I  detrend  the  population  patterns  by  a  linear 
regression  against  the  spatial  gradient.  Let  px  be  the  vector  of  prey  or  predator  values  in 
space  a;  at  a  fixed  time.  I  assume  that  px  is  the  following  linear  function  of  the  gradient  rx, 


Px  —  co  T  fdfx  -f-  €x , 


(5.4) 


where  the  residuals  ex  are  uncorrelated  random  errors  with  mean  0  and  equal  variances  in 
space.  The  parameters  a  and  f3  are  estimated  by  least  squares.  The  characteristic  scale 
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(i.e.  the  correlation  length)  is  then  computed  for  the  residuals.  This  allows  comparisons 
of  characteristic  scales  among  linearly  detrended  population  patterns  for  different  diffusion 
coefficients. 

Results 

Figure  5-2  shows  the  distributions  of  cross-correlation  coefficients  between  predator  patterns 
and  the  underlying  gradient  for  different  values  of  the  diffusion  coefficient  (i.e.  for  hc,  hq 
and  hp).  For  the  larger  value  of  diffusion,  when  the  dynamics  is  periodic,  the  predator 
patterns  strongly  reflect  the  gradient  (Figure  5-2,  C).  In  the  route  to  chaos,  as  the  diffusion 
coefficient  decreases,  predator  patterns  differ  more  and  more  from  the  environmental  pattern 
(Figure  5-2,  A,B).  A  similar  result  holds  for  the  prey  (Figure  5-3).  For  the  larger  value  of 
diffusion,  the  prey  patterns  resemble  the  gradient  but  can  display  the  opposite  slope  (Figure 
5-3,  C).  For  low  diffusion  values,  when  the  dynamics  is  chaotic,  the  distribution  clusters 
around  zero,  indicating  that  frequently  in  time,  the  prey  displays  low  or  no  correlation  to 
the  environmental  gradient  (Figure  5-3,  A). 
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Figure  5-2:  Histograms  of  cross-correlation  coefficients  between  predator  patterns  and  spa¬ 
tial  gradient. A:  d  =  10~ 4 ,  chaotic  dynamics;  B:  d  =  10~3,  quasiperiodic  dynamics;  C: 
d  =  10-2,  periodic  dynamics. 
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Figure  5-3:  Histograms  of  cross- correlation  coefficients  between  prey  patterns  and  spatial 
gradient.  A:  d  =  10-4,  chaotic  dynamics;  B:  d  =  10~3,  quasiperiodic  dynamics;  C  :d  =  10~2, 
periodic  dynamics. 
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Figure  5-4:  Histograms  of  cross-correlation  coefficients  between  predator  patterns  and  spa¬ 
tial  gradients.  The  diffusion  coefficient  is  the  same  for  all  graphs.  A  and  B  differ  in  the 
slope  of  the  gradient  (in  A:  r(x)  =  —  1 .4m  +  2;  in  B:  r{x )  =  — 0.1m  +  0.7).  A  and  C  differ  in 
the  local  dynamics  that  diffusion  couples  (in  A:  predator  mortality  m  =  0.6  and  diffusion 
couples  local  limit  cycles;  in  B:  m  =  0.68  and  diffusion  couples  local  equilibria). 
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Is  this  low  correlation  value  just  a  result  of  the  low  diffusion  coefficient?  Figure  5-4 
shows  the  distributions  of  cross-correlation  coefficients  for  different  qualitative  dynamics 
but  same  diffusion  coefficient  (i.e.  for  hc,  hgg  and  he).  When  the  system  converges  to  a 
steady-state,  predator  pattern  perfectly  reflects  the  gradient  ;  while  for  chaotic  dynamics, 
there  is  frequently  low  correlation  to  the  environment  (Figure  5-4,  A,C).  The  steepness  of 
the  gradient  also  influences  the  cross-correlation  coefficient  (compare  Figures  5-4,  A  and 
B). 

The  histograms  of  characteristic  scales  provide  further  evidence  for  the  increasing  mist- 
match  between  population  and  environmental  patterns  as  the  diffusion  coefficient  decreases 
(Figures  5-5  and  5-6).  Note  that  in  these  plots,  I  have  given  the  arbitrary  characteristic 
scale  of  one,  to  any  pattern  whose  correlogram  failed  to  cross  zero  in  (0, 1).  This  allows  me 
to  include  in  the  distribution,  patterns  that  basically  reflect  the  gradient.  For  instance,  for 
the  higher  diffusion  coefficient  and  periodic  dynamics  (hp),  the  predator  always  displays  a 
trend  similar  to  the  gradient  (Figure  5-5,  C).  In  the  route  to  chaos,  as  the  diffusion  coeffi¬ 
cient  decreases,  long  trends  become  less  frequent  and  patterns  present  smaller  characteristic 
scales  (Figure  5-5,  A,B).  This  decrease  in  Lc  does  not  result  only  from  a  low  diffusion  co¬ 
efficient.  A  steeper  gradient  leads  to  shorter  scale  values  for  the  same  diffusion  coefficient 
(compare  Figures  5-6,  A  and  B).  Equilibrium  dynamics  for  a  low  diffusion  value,  produces 
a  long  trend  in  the  predator  pattern  (Figure  5-6,  C).  Finally,  the  histograms  of  residuals 
confirm  the  presence  of  shorter  spatial  scales  for  weak  diffusion  when  dynamics  are  chaotic 
(Figure  5-7). 

5.2.2  Complexity  of  population  patterns 

Solutions  were  investigated  so  far  by  snapshots  of  spatial  patterns  at  fixed  times.  Here,  I 
would  like  to  compare  the  complexity  of  the  whole  spatio-temporal  solution  for  different 
diffusion  coefficients.  One  possible  approach  to  quantify  complexity  is  to  estimate  the  at¬ 
tractor’s  dimension  (see  Chapter  4).  I  take  below  a  different  approach  that  is  more  easily 
related  to  the  spatio-temporal  properties  of  the  solutions.  Results  will  confirm  the  appear¬ 
ance  of  smaller  spatial  scales  and  more  elaborate  patterns  as  diffusion  becomes  weaker. 
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Figure  5-5:  Histograms  of  characteristic  scales  for  predator  patterns.  A:  d  =  10  4,  chaotic 
dynamics;  B:  d  =  10~3,  quasiperiodic  dynamics;  C :  d  =  10~2,  periodic  dynamics. 
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Figure  5-6:  Histograms  of  characteristic  scales  for  predator  patterns.  The  diffusion  co¬ 
efficient  is  the  same  for  all  graphs.  A  and  B  differ  in  the  slope  of  the  gradient  (in  A: 
r(x)  =  -1.4a:  +  2;  in  B:  r(x)  =  -0.1a:  +  0.7).  A  and  C  differ  in  the  local  dynamics  that 
diffusion  couples  (in  A:  predator  mortality  m  —  0.6  and  diffusion  couples  local  limit  cycles; 
in  B:  m  =  0.68  and  diffusion  couples  local  equilibria). 
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Figure  5-7:  Histograms  of  characteristic  scales  for  predator  patterns  with  linear  trend  re¬ 
moved.  A:  d  =  10~4,  chaotic  dynamics;  B:  d  =  10~3,  quasiperiodic  dynamics;  C:  d  =  10~2, 
periodic  dynamics. 

Method 

Let  y(t,  x )  be  the  M  X  N  matrix  of  population  values  with  rows  corresponding  to  time  and 
columns  corresponding  to  space.  The  goal  is  to  approximate  y(x,t )  as  the  following  sum, 

y(x,  t)  «  Bn(t)Vn(x),  (5.5) 

n 

where  the  functions  Vn  depend  only  on  space  and  the  coefficients  B„  only  on  time.  I 
propose  to  measure  complexity  by  the  number  of  spatial  modes  Vn  needed  to  approximate 
the  solution  to  a  given  degree  of  accuracy. 

I  obtain  the  sum  of  equation  5.5  by  a  method  familiar  to  ecologists  as  Principal  Com¬ 
ponent  Analysis  (PCA)  and  to  oceanographers  as  Empirical  Orthogonal  Functions  (EOF). 
This  method  is  currently  used  in  turbulence  to  identify  coherent  structures  in  turbulent 
flows  and  to  study  the  attractors  of  spatio-temporal  dynamical  systems  (Sirovich,  1987; 
Sirovich  and  Rodriguez,  1987).  It  has  been  recently  extended  to  simplify  models  that  cou¬ 
ple  biology  and  physics  in  biological  oceanography  (Flierl  and  Davis,  in  press).  It  is  based 
on  the  orthogonal  decomposition  of  a  covariance  matrix. 
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The  spatial  functions  Vn  are  obtained  as  the  eigenvectors  of  the  N  X  N  covariance  matrix 


M  =  yTy  (5.6) 

where  yT  denotes  the  transpose  of  y.  Note  that  M  is  a  symmetric  matrix  and  therefore 
the  eigenvectors  Vn  are  real  and  orthogonal.  I  let  V\  denote  the  eigenvector  with  the 
dominant  eigenvalue,  V2,  the  eigenvector  with  the  subdominant  eigenvalue,  and  so  on.  The 
coefficients  Bn  are  obtained  by  projecting  the  solution  onto  the  new  coordinate  system 
provided  by  the  eigenvectors  Vn.  First  select  a  number  of  modes  Vn  and  form  the  matrix 
wn  =  [Vi  V2....V„].  Then,  the  coefficients  Bn  correspond  to  the  columns  of  the  matrix 
yW%  and  the  approximation  in  equation  5.5  is  given  by  BWn. 

I  evaluate  the  accuracy  of  the  approximation  with  a  mean  square  difference  that  includes 
all  temporal  and  spatial  locations.  Because  the  functions  Vn  condense  the  spatial  patterns 
for  the  whole  time  span  considered,  the  number  n  required  for  a  good  approximation  of  the 
solution  provides  a  measure  of  spatio-temporal  complexity. 

Results 

I  compare  the  predator  spatio-temporal  patterns  at  three  different  diffusion  coefficients  (i.e. 
hc,  hq,  and  hp).  For  weak  diffusion  (hc),  seven  modes  are  needed  for  a  high  degree  of 
accuracy  (mean  square  error=6  X  10~4)  (Figure  5-8).  A  subset  of  the  spatial  functions  Vn  is 
shown  in  Figure  5-9.  The  function  Vj  captures  the  spatial  gradient  and  successive  modes  Vn 
capture  patterns  with  increasingly  smaller  scales.  For  a  larger  diffusion  coefficient  ( hq ),  three 
modes  are  needed  for  a  similar  degree  of  accuracy  (mean  square  error  =  4.5  X  10~4)  (Figure 
5-10).  For  even  stronger  diffusion  (hp),  two  modes  approximate  the  solution  extremely  well 
(mean  square  error  =  2  X  10-5)  (Figure  5-11).  For  both  cases  ( hq  and  hp),  the  functions  Vf 
(not  shown  here)  represent  the  long  linear  trend  in  the  solution.  Fewer  additional  functions 
are  needed  as  the  spatial  patterns  more  closely  resemble  the  gradient.  Thus,  in  the  route  to 
chaos,  as  the  diffusion  coefficient  decreases,  spatio-temporal  solutions  become  more  complex 
and  display  smaller  spatial  scales. 

Up  to  this  point,  I  have  focused  on  comparing  the  pattern  of  the  linear  environmental 
gradient  to  those  of  the  predator  and  prey.  I  next  investigate  other  spatial  properties  of  the 
gradient  and  how  they  affect  the  spatio-temporal  dynamics  of  the  system. 
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Figure  5-8:  Approximation  to  chaotic  solutions.  Temporal  solutions  for  predator  numbers 
are  shown  at  a  fixed  location  in  space,  x  =  0.85,  for  d  =  10-4  (continuous  line).  Approx¬ 
imations  (dotted  line)  use  an  increasing  number  of  modes  from  top  to  bottom  (A:  three 
modes;  B:  five  modes;  C:  seven  modes). 
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Figure  5-9:  Orthogonal  spatial  functions  for  chaotic  predator  solutions  (d  =  10-4).  Only 
four  of  the  seven  functions  used  in  the  approximation  of  figure  5-8  are  shown.  See  text  for 
details. 


118 


Time 


Time 


Figure  5-10:  Approximation  to  quasiperiodic  solutions.  Temporal  solutions  for  predator 
numbers  are  shown  at  a  fixed  location  in  space,  x  =  0.85,  for  d  =  10~3  (continuous  line). 
Approximations  (dotted  line)  use  an  increasing  number  of  modes  from  top  to  bottom  (A: 
one  mode;  B:  two  modes;  C:  three  modes). 
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Figure  5-11:  Approximation  to  periodic  solutions.  Temporal  solutions  for  predator  numbers 
are  shown  at  a  fixed  location  in  space,  x  =  0.85,  for  d  =  10-2  (continuous  line).  Approxi¬ 
mations  (dotted  line)  use  an  increasing  number  of  modes  from  top  to  bottom  (A:  one  mode; 
B:  two  modes). 
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5.3  On  essential  properties  of  the  spatial  gradient 

5.3.1  Limit  cycles 

Recall  that  the  parameters  of  the  linear  gradient  were  chosen  to  obtain  limit  cycles  of  the 
predator  and  prey  at  each  fixed  spatial  location  in  the  absence  of  diffusion.  Are  these 
local  limit  cycles  necessary  for  complex  dynamics?  It  is  well  known  that  the  interaction  of 
two  or  more  oscillators  in  time  can  lead  to  complex  dynamics  via  the  quasiperiodic  route 
to  chaos  (for  ecological  examples  see  Kot  et  al .,  1992  ;  Hastings  and  Powell,  1991).  By 
analogy,  complex  dynamics  in  system  4.3  would  result  from  the  spatial  coupling  of  local 
predator-prey  cycles.  In  fact,  other  examples  of  diffusion-induced  chaos  have  been  explained 
as  the  result  of  the  spatial  coupling  of  local  oscillators  (Sharrett  et  al .,  1995;  Vastano  et 
al. ,  1990).  In  the  model  studied  here,  the  local  limit  cycles  differ  in  frequency  because  of 
the  underlying  gradient  (Figure  5-12,  curve  for  d  =  0).  To  examine  the  frequency  structure 
resulting  from  weak  diffusion,  I  computed  a  mean  period  of  oscillation  at  a  fixed  location  as 
the  mean  time  between  successive  local  maxima.  Figure  5-12  reveals  the  existence  of  spatial 
intervals  in  which  the  mean  period  of  oscillation  remains  constant.  These  frequency-locked 
regions  could  act  as  local  oscillators  whose  interplay  drives  the  chaotic  dynamics.  Note 
that  oscillations  with  a  rational  frequency  would  yield  a  return  map  of  local  maxima  that  is 
periodic,  and  therefore  given  by  a  finite  set  of  points.  By  contrast,  the  mean  frequency,  at 
any  location  within  a  region  in  figure  5-12,  appears  irrational:  the  return  map  is  aperiodic 
and  points  of  the  map  are  never  revisited.  As  a  consequence,  the  frequencies  of  two  adjacent 
regions  would  be  incommensurate.  It  is  known  that  in  temporal  systems,  the  interplay  of 
two  oscillators  with  incommensurate  frequencies  leads  to  complex  dynamics  (Kot  et  al, 
1992).  (As  a  brief  speculation,  the  step  pattern  of  this  curve  is  reminiscent  of  the  Cantor 
function  known  as  the  devils’s  staircase.) 

The  importance  of  the  local  limit  cycles  is  further  supported  by  a  large  number  of 
simulations  that  failed  to  produce  complex  dynamics  when  diffusion  coupled  local  equilibria. 
For  example,  for  a  higher  predator  mortality  rate  and  weak  diffusion,  the  system  converges 
to  a  steady-state  (see  solutions  he  and  pe  in  section  5.2.1). 

Diffusion,  however,  can  lead  to  chaos  when  limit  cycles  occur  only  in  part  of  the  gradient. 
To  show  this,  I  chose  a  gradient  that  produces  both  limit  cycles  and  equilibria  in  the 
absence  of  diffusion.  I  incorporated  the  spatial  gradient  in  the  parameter  a,  instead  of  r, 
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Figure  5-12:  Mean  period  of  oscillation  as  a  function  of  space.  Without  diffusion,  the  local 
dynamics  is  periodic.  The  period,  which  varies  in  space  because  of  the  gradient,  is  plotted  as 
a  function  of  space  (continuous  line,  d  =  0).  For  weak  diffusion,  when  the  system  is  chaotic, 
the  mean  period  of  oscillation  at  a  given  location  is  obtained  as  the  mean  time  between 
successive  local  maxima.  Its  distribution  in  space  displays  a  staircase  pattern  (dotted  line, 
d  =  10~4).  The  return  map  for  the  local  maxima  is  shown  for  two  different  locations. 
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the  non-dimensional  growth  rate  of  the  prey.  Recall  that  the  non-dimensional  parameter 
a  =  ( CiK)/(C2R )  where  C\,  C2 ,  K  and  R  are  dimensional  parameters  of  the  model  4.3 
denoting  respectively,  the  predator’s  maximum  uptake  rate  and  half-saturation  constant  of 
uptake,  and  the  prey’s  carrying  capacity  and  maximum  growth  rate.  Thus,  a  gradient  in 
the  maximum  uptake  rate  of  the  predator  would  generate  a  gradient  in  a.  This  choice  of 
gradient  is  necessary  because  the  value  of  a  determines  the  local  stability  of  the  equilibrium 
for  the  predator-prey  interaction  in  equations  4.3,  while  the  value  of  r  does  not.  I  therefore 
considered  the  predator-prey  system 


dp 

dt 

dh 

dt 


rp(  1  -  P) 
axp 


axP  u  ,  jd2P 
1  +  bp  +  dx2' 


,  .  ,d2h 

l  +  bph-mh  +  d8-I i- 


(5.7) 


where  the  nondimensional  parameters  are  defined  as  before,  but  r  —  1  and  ax  =  e  +  fx.  In 
the  absence  of  diffusion,  at  a  fixed  point  in  space,  the  predator-prey  interaction  of  system 
5.7  has  a  locally  stable  positive  equilibrium  when 


a  < 


bm{b  +  1) 
6-1 


a  >  6m,  and  6  >  1, 


(5.8) 


(see  Appendix). 

In  the  absence  of  diffusion,  for  parameters  e  =  9 ,  /  =  -6,  6  =  3,  and  m  =  0.6,  the 
predator-prey  system  displays  limit  cycles  in  x  =  [0,0.9)  and  equilibria  in  x  =  [0.9,1].  In 
spite  of  the  local  equilibrium  dynamics,  diffusion  can  induce  chaos  in  this  system.  Figure 
5-13  illustrates  the  irregular  spatio-temporal  patterns  of  the  prey  in  the  chaotic  regime 
(d  =  10~4).  The  temporal  population  patterns  appear  highly  irregular  but  are  strongly 
damped  in  amplitude  at  the  end  of  the  gradient  where  diffusion  couples  equilibria  (Figure 
5-14).  The  chaotic  nature  of  the  dynamics  is  demonstrated  by  its  sensitivity  to  initial 
conditions  (i.e.  a  positive  Lyapunov  exponent  A  =  1.5  X  10-2  bits  per  unit  time). 

5.3.2  Steepness  of  the  gradient 

The  steepness  of  the  gradient  is  an  important  determinant  of  the  dynamics.  This  is  not 
surprising  since  the  slope  of  the  gradient  determines  the  range  of  frequencies  and  amplitudes 
of  the  predator-prey  cycles  that  diffusion  couples. 
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Figure  5-13:  Spatio-temporal  prey  patterns.  Diffusion  couples  limit  cycles  in  the  interval 
x  —  [0,0-9)  and  equilibria  in  x  =  [0.9, 1). 

Figure  5-15  demonstrates  the  changes  in  dynamics  for  different  slopes  of  the  gradient 
rx  but  same  diffusion  coefficient  d  =  10-4.  This  bifurcation  diagram  is  obtained  as  fol¬ 
lows:  a  spatial  location  is  chosen,  and  then,  the  successive  local  maxima  in  prey  numbers 
at  this  fixed  point  are  plotted  for  the  corresponding  slope  value.  A  period  1  trajectory  pro¬ 
duces  a  single  point.  More  generally,  periodic  trajectories  produce  finite  number  of  points. 
Successive  maxima  of  quasiperiodic  and  chaotic  trajectories  spread  over  a  range  of  values. 
Whereas  the  former  densely  cover  this  range,  the  latter  exhibit  more  structure.  For  a  de¬ 
tailed  description  of  the  changes  in  dynamics  associated  with  different  qualitative  regions  in 
a  bifurcation  diagram  of  the  model,  see  Chapter  4.  Figure  5-15  shows  that,  as  the  absolute 
value  of  the  slope  increases,  the  model  dynamics  change  from  periodic  to  aperiodic.  For  a 
slope  equal  to  —1.4,  the  parameters  correspond  to  the  study  in  Chapter  4,  where  I  have 
shown  that  the  dynamics  are  chaotic. 


5.3.3  Nonlinear  gradients 

This  work  has  so  far  focused  on  a  linear  gradient.  Figure  5-16  shows  that  complex  dy¬ 
namics  also  occurs  for  nonlinear  gradients.  The  Poincare  sections  in  the  top  panels  are 
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Figure  5-14:  Temporal  behavior  of  predator  numbers.  Two  spatial  locations  in  figure  5-13 
are  chosen  to  illustrate  the  long-term  temporal  dynamics  of  the  system  (A:  x  =  0.95;  B: 
x  =  0.15).  Predator  patterns  are  irregular  at  both  locations  but  appear  strongly  damped 
in  amplitude  at  the  end  of  the  gradient  where  diffusion  couples  equilibria. 

obtained  by  first  reconstructing  the  attractor  in  three  dimensions  from  the  time  series  of 
prey  numbers  at  a  fixed  location,  and  then  cutting  the  reconstructed  trajectory  with  a  plane 
(see  Chapter  4).  For  periodic  dynamics,  the  Poincare  section  is  a  finite  set  of  points;  for 
quasiperiodic  dynamics,  a  one-dimensional  curve  topologically  equivalent  to  a  circle.  The 
higher-dimensional  attractor  in  Figure  5-16  (top  panels)  is  indicative  of  chaos.  Further 
evidence  is  provided  by  sensitivity  to  initial  conditions  (i.e.  positive  Lyapunov  exponents, 
A  =  5.7  x  10~3  (A)  and  A  =  1.8  x  10-3  bits  per  unit  time  (B)). 

These  examples  consider  gradients  given  by  long  spatial  trends.  To  examine  the  effect 
of  spatial  heterogeneity  with  higher  frequency  variation,  I  considered  gradients  of  the  form 

0.5Af(3  -  cos(2ttN(x  —  1))),  (5-9) 


that  is,  sinusoidal  curves  with  maxima  2 M,  minima  M,  and  N  peaks  in  x  -  (0, 1).  Examples 
of  gradients  with  the  same  M  but  different  N  are  shown  in  Figure  5-17  (A).  Notice  that  as 
N  increases,  the  spatial  frequency  of  the  gradient  increases.  This  effectively  increases  the 
average  absolute  value  of  the  local  gradient,  and  decreases  the  distance  between  any  two 
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Figure  5-15:  Bifurcation  diagram  for  an  increasing  slope  of  the  gradient.  Prey  maxima  at  a 
fixed  location  in  space  (x  =  0.85)  are  plotted  for  different  values  of  the  slope  of  the  gradient 


frequencies  in  the  predator-prey  cycles  that  are  coupled  by  diffusion.  These  effects  are  also 
produced  by  increasing  the  slope  of  a  linear  gradient.  For  the  sinusoidal  gradients,  however, 
the  overall  range  of  local  frequencies  in  the  predator-prey  interaction  is  not  modified.  Figure 
5-17  demonstrates  that  complex  dynamics  requires  a  weaker  and  weaker  diffusion  coupling 
as  the  frequency  of  the  gradient  increases.  I  plot,  for  different  N,  the  diffusion  coefficient 
d  at  which  the  system  changes  qualitative  dynamics  from  periodic  to  quasiperiodic.  This 
critical  d  value  shows  a  decreasing  trend  with  N. 


Figure  5-16:  Chaotic  dynamics  for  two  nonlinear  gradients.  Poincare  sections  obtained  from 
reconstructed  attractor  (top  panels)  and  the  corresponding  gradients  in  prey  growth  rate 
(bottom  panels). 
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Figure  5-17:  Bifurcation  to  complex  dynamics  for  nonlinear  gradientsSinusoidal  gradients 
are  shown  in  A  (diamonds:  N  =  2;  dotted  line:  N  =  1;  continuous  line:  N  =  0.5).  In  B, 
the  critical  value  of  d  at  which  bifurcation  from  periodic  to  quasiperiodic  dynamics  occur, 
is  plotted  for  sinusoidal  gradients  of  increasing  frequency  (i.e.,  increasing  N). 
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5.3.4  Diffusion  rates  of  predator  and  prey 

All  simulations  in  this  study  considered  equal  diffusion  rates  of  predator  and  prey.  However, 
for  a  given  environmental  gradient  predator  and  prey  may  diffuse  at  different  rates.  Are 
equal  diffusion  rates  of  predator  and  prey  a  necessary  requirement  for  complex  dynamics? 
To  address  this  question,  I  varied  the  diffusion  rate  of  the  predator  for  the  same  diffusion 
rate  of  the  prey  (d  =  10~4).  The  resulting  bifurcation  diagram  is  shown  in  Figure  5-18. 
This  figure  resembles  qualitatively  the  bifurcation  diagram  obtained  in  Chapter  4  by  varying 
both  diffusion  coefficients  when  they  are  equal.  As  d  decreases,  the  model  dynamics  changes 
from  periodic  to  quasiperiodic,  to  chaotic  dynamics.  A  requirement  for  complex  dynamics 
is  therefore  that  both  species  diffuse  at  low  but  not  necessarily  equal  rates. 


0.20  1 - - - 1 - 1 - 1 - * - 1 - 1 - 1 
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Figure  5-18:  Bifurcation  diagram  for  a  decreasing  diffusion  coefficient  of  the  predator.  Prey 
maxima  at  a  fixed  spatial  location  x  =  0.85  are  shown  as  a  function  of  the  predator’s 
diffusion  coefficient.  The  prey’s  diffusion  coefficient  remains  the  same  ( d  =  10~4). 


129 


5.4  Discussion 


In  the  first  part  of  this  chapter,  I  have  shown  that  the  spatial  patterns  of  a  predator  and  its 
prey  diffusing  in  heterogeneous  space,  can  differ  strongly  from  the  underlying  environmental 
gradient.  In  the  route  to  chaos,  as  diffusion  becomes  weaker,  this  difference  is  magnified 
and  the  predator  and  prey  patterns  display  smaller  spatial  scales. 

How  do  these  results  compare  to  the  recent  studies  by  McLaughlin  and  Roughgarden 
(1991,  1992)  of  predator-prey  interactions  and  diffusion  in  heterogeneous  space?  They 
considered  reaction- diffusion  equations  with  Lotka-Volterra  type  terms  for  the  predator-prey 
interaction,  and  focused  on  the  differential  mobility  of  predator  and  prey.  They  showed  that 
the  predator-prey  interaction  ‘sharpens’  the  underlying  environmental  pattern  by  reflecting 
it  in  the  prey  distribution  to  a  degree  proportional  to  the  mobility  of  predators  relative  to 
prey.  ‘Sharpens’  means  that  the  prey  patterns  can  display  steeper  local  slopes  than  the 
gradient.  In  fact,  when  the  predator  diffuses  on  a  sinusoidal  gradient  much  faster  than 
its  prey,  the  prey  cannot  survive  in  regions  of  low  growth  rate  and  displays  steep  patterns 
in  patches  of  high  growth  rate  (McLaughin  and  Roughgarden,  1992).  By  contrast  to  the 
results  presented  here,  however,  the  population  patterns  have  the  same  dominant  spatial 
scale  than  that  of  the  underlying  gradient:  the  positions  of  local  maxima  in  the  gradient 
and  prey  numbers  always  coincide.  I  infer  that  the  correlation  coefficient  between  gradient 
and  population  patterns  would  generally  be  high  in  their  model,  except  for  predators  that 
diffuse  extremely  fast  relatively  to  the  prey.  These  differences  with  the  patterns  obtained 
here  may  result  from  the  different  diffusion  ranges  considered:  in  their  work,  predators 
diffuse  faster  than  the  prey;  here,  diffusion  was  weak  for  both  predator  and  prey.  Another 
possible  explanation  regards  the  different  formulations  of  the  predator-prey  interaction:  the 
Lotka-Volterra  equations  display  neutral  cycles  (i.e.  cycles  that  are  structurally  unstable 
since  any  perturbation  leads  to  a  cycle  with  different  amplitude);  the  predator-prey  model 
used  here  diplays  limit  cycles.  It  is  an  open  question  whether  the  diffusive  coupling  of  these 
neutral  cycles  can  lead  to  complex  dynamics.  In  their  simulations,  the  coupling  of  neutral 
cycles  leads  to  periodic  or  equilibrium  solutions. 

In  a  second  part  of  this  chapter,  I  have  explored  properties  of  the  gradient  that  influence 
the  spatio-temporal  dynamics  of  the  predator  and  prey.  Results  suggest  that  both  the 
gradient  and  the  local  limit  cycles  are  required  for  complex  dynamics  in  the  model.  Complex 
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dynamics  would  result  from  the  spatial  coupling  by  diffusion  of  local  oscillators  that  differ 
in  frequency  because  of  the  underlying  gradient.  The  importance  of  local  limit  cycles  is 
supported  by  studies  of  diffusion-driven  instabilities  and  chaos  in  physical  and  chemical 
systems  (Nicolis  and  Gaspard,  1990;  Vastano  et  al.,  1990).  Recent  work  by  Sherratt  et 
al.  (1995)  indicates,  however,  that  environmental  heterogeneity  is  not  a  requirement  for 
spatio-temporal  chaos  in  reaction-diffusion  models  of  predator-prey  interactions.  In  their 
model,  chaos  is  generated  in  the  wake  of  invasive  waves  of  predators.  A  wave  of  this  type 
develops  when  predators  are  introduced  locally  into  a  uniform  distribution  of  the  prey.  This 
explains  why  this  behavior  was  not  observed  here  with  the  more  general  initial  conditions 
of  predator  and  prey  present  in  the  whole  domain.  It  is  possible  that  the  predator  wave 
in  their  model  effectively  generates  a  local  gradient.  Two  other  ecological  examples  of 
chaos  in  reaction-diffusion  systems  with  homogeneous  environments  regard  a  predator  and 
two  competing  prey  and  a  network  of  predator-prey  interactions  (Ikeda  and  Mimura,  1993; 
Pahl-Wostl,  1993). 

Spatial  models  will  provide  useful  tools  to  investigate  how  methods  to  detect  nonlinear¬ 
ity,  dimensionality,  and  sensitivity  to  initial  conditions,  perform  under  different  sampling 
regimes.  In  ecology,  such  methods  are  usually  tested  with  low-dimensional  temporal  sys¬ 
tems.  Spatial  interactions  introduce  a  higher  dimensionality.  Because  ecological  time  series 
are  generally  short,  estimates  of  Lyapunov  exponents  would  be  local  in  both  space  and 
time  (for  a  definition  of  local  Lyapunov  exponents  see  Ellner  and  Turchin,  1995).  In  addi¬ 
tion,  ecological  data  sets,  particularly  oceanographic  ones,  reflect  the  mixing  of  space  and 
time  introduced  by  the  sampling  scheme.  Little  at  al.  (in  preparation)  have  begun  to  ex¬ 
plore  the  consequences  of  Lagrangian  vs.  Eulerian  sampling  for  the  detection  of  chaos  and 
nonlinearity.  This  issue  is  specially  relevant  for  plankton  systems. 

The  qualitative  results  of  this  chapter  would  apply  to  planktonic  organisms  transported 
by  turbulent  diffusion  (Pascual,  1993).  It  would  be  interesting,  however,  to  explore  diffusion- 
induced  complex  dynamics  in  a  model  specifically  tailored  to  a  planktonic  system.  In¬ 
deed,  Malchow  (1993)  comments  on  the  occurrence  of  apparently  complex  dynamics  in  a 
phytoplankton-zooplankton  reaction-diffusion  model  that  incorporates  a  time-  and  depth- 
dependent  growth  rate  of  the  prey.  The  resulting  spatial  patterns  continuously  change.  I 
suggest  that  those  patterns  would  differ  in  spatial  scale  from  the  underlying  gradient  in 
light,  and  therefore,  that  a  linear  approach  relying  on  scale  matching  would  fail  to  reveal 
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their  cause. 


Although  biological  oceanography  has  pioneered  the  concept  of  scales  in  ecology  (Haury 
et.,  1978;  Steele,  1978),  a  linear  perspective  has  dominated  the  study  of  how  scales  integrate 
(Denman,  1994;  Denman  and  Powell,  1984).  A  few  authors  have  cautioned  against  simple 
assumptions  of  linearity  (Dwyer  and  Perez,  1983;  Denman  and  Powell,  1984;  Star  and 
Cullen,  1981,  Steele,  1988).  This  work  further  supports  the  need  for  a  nonlinear  perspective. 


5.5  Appendix:  Local  stability  of  positive  equilibrium 


Consider  the  predator  system  (5.7)  with  no  diffusion  and  no  spatial  gradient  (i.e.  ax  =  a). 
The  nontrivial  equilibrium  is  given  by 


{.Pei  he) 


(  rn  r(l  -  pe)(l  +  bpe) 
\a  —  bm  ’  a 


(5.10) 


which  is  positive  provided  a  >  bm.  Linearization  at  this  equilibrium  gives  the  following 
Jacobian  matrix, 

j  =  (  r<l  -  2f«)  “  tUvF  -»).  (5.11) 

\  (ufep  0  / 

Conditions  for  the  local  stability  of  the  equilibrium  are: 

1.  -Trace(J)  >  0 


2.  Determinant /)  >  0 

(May,  1973).  Condition  (2)  is  always  satisfied.  Condition  (1)  is  equivalent  to 

a(b  -  1)  <  bm{b  +  1).  (5.12) 

Notice  that  this  inequality  imposes  no  condition  on  r,  but  does  on  a  for  given  b  and  m. 
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Chapter  6 


Intermittency  in  the  plankton:  a 
multifractal  analysis  of 
zooplankton  biomass  variability- 


...j’ai  dit  souvent  que  les  branches  des  arbres  etaient 
elles-memes  de  petits  arbres  complets:  des  fragments 
de  rochers  sont  semblables  a  des  masses  de  rochers, 
des  particules  de  terre  a  des  amas  enormes  de  terre. 

Je  suis  persuade  qu’on  trouverait  en  quantite  de  ces 
analagies.  Une  plume  est  composee  d’un  million  de  plumes... 

— Delacroix.  Journal,  1823-66 


6.1  Introduction 

Plankton  data  vary  on  a  wide  range  of  spatial  and  temporal  scales  (Haury  et  al.,  1978;  Steele, 
1978).  Describing  this  variability  is  an  important  problem  in  plankton  ecology,  especially 
given  recent  developments  in  methods  for  continuously  recording  data  at  high  spatial  and 
temporal  resolution  (Dickey,  1988,  1991).  Quantitative  characterization  of  pattern  provides 
a  basis  for  comparing  models  to  data,  and  biological  to  environmental  fluctuations.  A  well 
known  approach  to  such  characterization,  spectral  analysis,  was  pioneered  in  ecology  by 
plankton  studies  (Platt  and  Denman,  1975). 

In  this  paper  we  explore  an  alternative  approach,  characterizing  zooplankton  biomass 
variability  as  a  multifractal.  Multifractals  are  a  generalization  of  fractals  (e.g.,  Mandlebrot 

1This  chapter  is  in  press  in  Pascual,  M.,  F.A.  Ascioti,  and  H.  Caswell  (1995)  J.  of  Plankton  Res.  17(5). 
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1983)  from  the  description  of  geometrical  patterns  to  the  description  of  spatial  or  temporal 
series  of  numerical  quantities.  The  basic  idea  of  fractal  pattern  is  that  a  power  law  describes 
the  relation  between  some  quantity  and  the  scale  on  which  it  is  measured.  The  exponent 
of  the  power  law,  known  as  the  fractal  dimension,  shows  how  the  quantity  relates  to  the 
scale  of  measurement.  A  well-known  example  is  the  problem  of  measuring  the  length  of  a 
coastline.  The  finer  the  scale  of  measurement,  the  longer  the  coastline  will  appear;  “the 
length”  of  the  coastline  is  not  a  well-defined  concept  (Richardson,  1961;  Mandelbrot,  1983). 
However,  the  variation  of  length  with  the  scale  of  measurement  is  well-described  by  a  power 
function.  This  provides  a  fractal  dimension  and  completely  characterizes  the  way  in  which 
the  length  of  the  coastline  varies  with  scale. 

Multifractals,  which  will  be  reviewed  below,  describe  patterns  by  scaling  relations  that 
require  a  family  of  different  exponents,  rather  than  the  single  exponent  of  fractal  patterns. 
They  are  particularly  well  suited  to  describing  quantities  that  vary  intermittently  (i.e., 
occasional  and  unpredictable  large  peaks  separated  by  very  low  values),  and  have  been 
applied  to  a  variety  of  intermittent  measures  associated  with  nonlinear  phenomena  in  physics 
and  geophysics  (Meakin,  1983;  Meneveau  and  Sreenivasan,  1991;  Prasad  et  al.  1988,  Ladoy 
et  a/.,  1991;  Sreenivasan,  1991). 

We  will  present  evidence  here  for  the  multifractal  structure  of  zooplankton  biomass 
variability.  Our  analysis  is  based  on  an  hourly  time  series  of  vertically  integrated  acoustic 
biomass  measurements,  taken  from  a  fixed  mooring  on  the  Atlantic  coastline.  We  analyzed 
two  estimates  of  variability:  the  first  difference  squared  and  the  squared  difference  from 
the  mean.  When  summed  over  time  these  quantities  provide  estimates  of  biomass  vari¬ 
ance.  Our  goal  is  to  describe  the  distribution  in  time  of  the  total  variability  in  the  data. 
This  distribution  is  highly  intermittent:  extreme  localized  contributions  account  for  a  large 
proportion  of  the  total  variability. 

6.2  The  data 

The  data  on  plankton  biomass  were  provided  by  C.  Flagg  of  the  Brookhaven  National  Lab¬ 
oratory.  Zooplankton  biomass  was  estimated  from  measurements  of  acoustic  backscatter 
intensity  at  a  fixed  mooring  off  the  continental  slope  of  Maryland.  Three  different  deploy¬ 
ments  of  an  Acoustic  Doppler  Current  Profiler  (ADCP)  operating  at  307.5  kHz  provided 
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three  time  series,  labeled  respectively  A,  B,  and  C.  The  instrument  recorded  data  from  10 
to  85  meters  of  depth,  at  3-minute  intervals  for  time  series  C  and  2.5-minute  intervals  for 
time  series  A  and  B  (Flagg  et  al. ,  1994).  Zooplankton  net  tows  were  used  to  calibrate  the 
instrument  and  convert  the  data  to  dry  weight  zooplankton  biomass  (mg/m3)  (Flagg  and 
Smith,  1989;  Flagg  et  al.,  1994).  The  data  analyzed  here  were  constructed  by  averaging 
measurements  vertically  and  hourly.  Figure  6-1  shows  the  resulting  time  series  of  zooplank¬ 
ton  biomass  obtained  at  the  three  different  deployments  from  5  February  1988  to  13  May 
1989.  The  series  contain  2732,  2735  and  4099  points,  respectively. 


TIME  (hours) 


Figure  6-1:  Time  series  of  zooplankton  biomass  dry  weight  (mg/m3).  (A):  Hourly  data 
from  February  16,  1988  (11:00  PM)  to  June  9,  1988  (6:00  PM).  (B):  Hourly  data  from  June 
28,  1988  (2:00  PM)  to  October  20,  1988  (12:00  AM).  (C):  Hourly  data  from  November  19, 
1988  (12:00  PM)  to  May  8,  1989  (6:00  PM). 


6.3  Methods 


6.3.1  The  multifractal  formalism 

Multifractal  analysis  requires  three  fundamental  terms:  support,  measure ,  and  measure 
density.  The  basic  data  consist  of  some  quantity,  which  we  will  refer  to  as  the  measure  (in 
our  case,  zooplankton  biomass  or  some  quantity  calculated  from  it)  along  an  axis  which 
we  will  refer  to  as  the  support  of  the  measure.  The  support  is  most  commonly  space  or 
time,  although  in  our  case  the  series  of  acoustic  measurements  contains  both  temporal  and 
spatial  components.  The  support  could  be  multidimensional,  although  in  our  case  it  is  a 
one- dimensional  axis.  The  measure  density,  as  its  name  suggests,  is  the  total  measure  over 
some  segment  of  the  support  axis,  divided  by  the  length  of  that  segment.  If  the  measure 
density  is  divided  by  the  total  measure  over  the  whole  data  set,  we  obtain  a  normalized 
measure  giving  the  proportion  of  the  total  measure  occuring  in  each  spatial  location. 

Multiplicative  processes 

To  introduce  the  concept  of  a  multifractal  measure,  we  consider  the  simplest  example  of 
a  process  generating  such  a  measure.  This  process,  known  as  the  self-similar  binomial 
process,  is  recursive.  One  starts  with  a  uniform  distribution  of  mass  over  the  unit  interval 
(0,1)  (Figure  2).  In  a  first  step  of  the  process,  the  unit  interval  is  subdivided  into  two  equal 
intervals  and  proportions  mo  and  1  —  mo  of  the  total  mass  (e.g.,  0.7  and  0.3  in  Figure  2) 
are  allocated  to  these  two  subintervals.  This  process  is  now  repeated  for  each  of  the  two 
subintervals:  they  are  subdivided  into  two  equal  intervals  and  their  corresponding  mass 
allocated  in  proportions  m0  and  1  -  m0  to  their  left  and  right  subintervals  respectively. 
Figure  6-2  shows  the  resulting  distribution  of  density  after  8  such  steps.  Notice  that  the 
measure  in  a  given  subinterval  is  the  integral  of  this  density. 

Two  important  properties  are  illustrated  by  this  example.  First,  the  density  is  intermit¬ 
tent,  infrequent  variations  of  large  amplitude  appear  within  more  regular  regions  of  lower 
values.  In  the  limit,  as  the  number  of  steps  in  the  binomial  process  becomes  arbitrarily 
large,  the  density  at  every  point  tends  to  either  zero  or  infinity.  At  the  points  where  the 
density  increases  without  bound  it  is  said  to  have  a  singularity.  Second,  the  measure  density 
exhibits  self-similarity  or  invariance  against  changes  of  scale,  in  the  following  sense.  After 
k  steps  of  the  process,  the  right  half  distribution  equals  the  left  half  times  m0/(  1  -  m0) 
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and  the  left  half  distribution  resembles  that  in  the  whole  interval  at  step  k  —  1  (see  Figure 
6-2).  In  fact  the  whole  distribution  can  be  obtained  from  the  left  half  by  stretching  it  in  the 
horizontal  direction  by  a  factor  of  2  and  in  the  vertical  direction  by  a  factor  of  1/(1  -  m0). 
As  the  numbers  of  steps  becomes  arbitrarily  large,  the  same  transformations  produce  the 
entire  distribution  from  its  left  half  portion.  Thus,  parts  of  the  distribution  resemble  the 
whole. 


Figure  6-2:  The  binomial  measure.  The  original  uniform  distribution  of  density  is  shown 
in  the  upper-left  panel.  The  first  8  fragmentation  steps  are  illustrated  in  following  panels 
(  counterclockwise  direction).  Notice  that  the  y-axis  corresponds  to  density  and  therefore, 
the  area  below  the  density  curve  provides  the  measure  in  any  given  subinterval. 
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The  binomial  process  is  a  special  case  of  a  larger  class  of  processes  called  multiplicative 
processes.  In  multiplicative  processes,  large  pieces  of  the  support  of  the  measure  break  down 
into  smaller  ones,  and  each  of  the  fragmented  pieces  yield  smaller  ones  and  so  on.  At  each 
step  of  this  cascade,  the  fragmented  pieces  receive  a  fraction  of  the  original  measure.  Thus, 
at  a  step  k  of  the  cascade,  the  measure  in  a  certain  fragment  will  be  given  by  the  product 
of  k  numbers  known  as  multipliers  (0.7  and  0.3.  in  the  example  above).  The  multipliers 
may  also  be  random  variables  with  a  certain  probability  distribution.  When  this  probability 
distribution  does  not  depend  on  the  step  of  the  cascade,  scale  similarity  results. 

Multiplicative  processes  provide  a  mathematical  ideal  of  multifractals  in  nature.  In  real 
data  sets  there  are  limits  to  the  scales  at  which  we  may  determine  a  measure,  to  the  number 
of  steps  that  a  multiplicative  cascade  can  achieve  and  to  the  range  of  scales  in  which  the 
multifractal  description  we  describe  below  would  apply. 


Describing  multifractal  processes 

An  intuitive  way  to  describe  a  measure  would  be  to  plot  the  frequency  distribution  of 
the  density;  i.e.,  a  distribution  showing  how  much  of  the  support  is  characterized  by  any 
specific  density.  However,  like  the  length  of  a  fractal  coastline,  the  frequency  distribution 
of  a  multifractal  density  changes  as  a  function  of  the  scale  of  measurement.  Therefore,  our 
attention  focuses  on  the  scaling  properties  of  the  measure. 

These  scaling  properties  require  not  one,  but  a  whole  family  of  different  exponents.  We 
present  below  two  families  of  exponents  of  which  the  first  has  inspired  the  name  multifrac¬ 
tals,  while  the  second  is  more  useful  for  analysis  of  empirical  data. 

Divide  the  support  of  a  measure  M  into  segments  of  length  r.  Let  Mr(x)  denote  the 
measure  in  one  such  fragment  centered  at  coordinates  x.  The  corresponding  density  is 
denoted  by  rar(x)  and  equals  Mr(x.)/rd  where  d  is  the  dimension  of  the  support  (i.e.  d  =  1 
for  a  time  axis  or  a  spatial  transect,  d  =  2  for  a  spatial  area,  etc...).  We  define  for  each 
segment  a  quantity  a(x)  defined  by 


a(x) 


In  {Mr /Ml) 
In  (r/L) 


where  L  is  the  length  of  the  total  support.  In  the  limit  as  r  goes  to  zero,  a  measures  the 
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scaling  of  the  measure  or  the  density  with  length  of  the  segment: 


Ml 


(6.1) 


Here  we  use  the  symbol  ~  to  mean  that  the  left-hand  side  approaches  a  constant  times  the 
right-hand  side  in  the  limit  of  small  r.  The  quantity  Mt/Ml,  the  measure  normalized  by 
its  total  value,  varies  between  0  and  1,  and  gives  the  proportion  of  the  total  measure  in  a 
segment  of  size  r  centered  at  x. 

The  local  exponent  a  describes  how  the  measure  and  the  density  change  with  changes 
in  the  length  r  of  the  segment  (technically,  f*(x)  measures  the  singularity  strength  of  the 
density  at  x).  Equation  6.1  shows  that  the  measure  increases  as  r  increases.  The  smaller  the 
value  of  a,  the  faster  this  increase  will  be  at  the  smallest  scales  r.  Thus,  as  a  decreases,  more 
and  more  of  the  measure  is  contributed  by  smaller  and  smaller  scales.  This  is  illustrated  in 
Figure  6-3  by  comparing  the  behavior  of  a  hypothetical  measure  at  two  different  points,  xi 
and  X2,  on  a  one  dimensional  support.  At  xt,  a(xi)  <  1  and  therefore,  Mr(xi)  increases 
rapidly  near  r  =  0;  at  X2,  a(x2)  >  1,  and  Mr(x 2)  varies  slowly  near  r  =  0  (Figure  6-3, 
b).  Correspondingly,  the  measure  displays  a  peak  and  a  trough  in  segments  of  the  support 
centered  at  xi  and  X2  (Figure  6-3,  a).  More  generally,  the  relation  between  a  and  the 
support  dimension  d  distinguishes  locations  of  the  support  with  high  ( a  <  d )  and  low 
(a  >  d)  local  intensity  of  the  measure.  These  correspond  to  locations  where  the  density 
grows  without  bounds  (a  <  d)  as  r  — *  0,  and  locations  where  the  density  approaches  zero 
( a  >  d)  as  r  decreases  (Equation  6.2).  The  smaller  the  value  of  a(x)  <  d,  the  sharper  the 
peak  in  the  density  at  location  x. 

Each  segment  of  length  r  is  now  associated  with  a  value  of  cc(x)  describing  how  the 
measure  changes  with  scale  around  x.  Let  Nr{a)  denote  the  number  of  intervals  of  length  r 
with  a  values  in  the  interval  (a,  a+da).  We  complete  our  characterization  of  the  multifractal 
by  seeing  how  Nr(a)  scales  with  r,  by  defining 


fr(a) 


log  Nr(a) 
log (r/L)  ' 


(6.3) 
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In  the  limit  as  r  goes  to  zero,  fr(a)  converges  to  a  well  defined  limit  /(a),  satisfying 

Nr(a)~(z)  (6.4) 

(Meneveau  and  Sreenivasan,  1991).  The  function  f(a)  is  called  the  singularity  spectrum. 
Because  of  the  similarity  of  expression  6.4  to  the  one  defining  fractal  dimension,  /(a) 
can  be  interpreted  as  the  fractal  dimension  of  the  set  of  intervals  with  a  in  (a,  a  +  da) 
(Frisch  and  Parisi,  1985).  Heuristically,  when  we  label  different  segments  of  size  r  with 
their  corresponding  a  value,  we  obtain  subsets  of  the  support  of  the  measure  made  of  all 
fragments  with  same  a.  These  subsets  are  geometrical  sets  and  in  the  limit,  as  r  becomes 
arbitrarily  small,  they  tend  to  sets  of  points.  Each  subset  has  its  corresponding  fractal 
dimension  /(a)  indicating  how  dense  it  is  in  the  measure  support.  If  f(a)  is  small,  the 
points  with  exponent  a  are  scattered  and  dustlike.  As  f(ot)  approaches  d,  the  set  of  points 
with  exponent  a  become  more  and  more  dense.  The  name  multifractal  results  from  having 
a  different  /(a)  for  each  a.  A  typical  parabolic  shape  of  the  singularity  spectrum  is  shown 
in  Figure  6-3(c). 

These  calculations  permit  us  to  define  precisely  the  notion  of  a  multifractal;  we  say  that 
a  pattern  is  multifractal  if  the  exponent  a,  defined  in  6.1,  spreads  over  a  range  of  values,  and 

r 

for  each  a  the  scaling  relation  6.4  holds.  The  variable  a  reflects  how  singular  (or  spiky)  the 
behavior  of  the  density  is  at  a  given  location,  and  its  corresponding  value  /(a)  how  frequent 
this  local  exponent  is  with  respect  to  other  values.  (For  a  more  detailed  discussion  of  when 
/(a)  can  be  interpreted  as  a  fractal  dimension  see  Meneveau  and  Sreenivasan,  1991).  The 
variation  in  a  values  is  characteristic  of  multifractal  measures;  exact  fractals,  by  contrast, 
have  the  same  a  for  all  locations  x.  The  variance  of  a  relates  to  the  degree  of  intermittency 
in  the  data  (see  section  6.5). 

Scaling  of  moments  in  multifractals 

A  second  way  to  characterize  multifractals,  and  one  which  is  readily  used  in  data  analysis, 
is  by  mean  of  moments.  Highly  intermittent  multifractal  signals  are  not  well  characterized 
by  a  few  low-order  moments,  such  as  the  ones  providing  the  mean  and  variance,  because  of 
the  strong  tails  of  their  probability  distributions.  Therefore,  the  approach  described  below 
relies  on  a  family  of  moments  and  their  respective  scaling  laws. 
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Figure  6-3:  The  exponents  a  and  f(a)  for  a  multifractal  measure  on  a  one  dimensional 
support.  In  (B),  the  behavior  of  the  measure  with  r  is  compared  at  two  different  points 
xi  and  x2.  At  xl5  a(xi)  <  1  and  therefore,  the  measure  Mr(x i)  increases  rapidly  near 
r  =  0;  at  x2,  a(x2)  >  1,  and  Mr(x2)  increases  slowly  near  r  =  0.  Correspondingly,  the 
measure  displays  local  high  and  low  intensity  in  segments  of  the  support  centered  at  xi  and 
x2  (A).  A  typical  parabolic  shape  of  the  singularity  spectrum  is  shown  in  (C).  When  f(a) 
is  small,  the  points  with  the  corresponding  exponent  a  are  scattered  and  dustlike.  As  f(a) 
approaches  d  =  1,  the  set  of  points  with  exponent  a  fill  more  and  more  of  the  support. 
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The  r/th  moment  of  a  quantity  x  can  be  denoted  by  ( xq ),  where  the  brackets  ()  denote 
the  expected  value.  For  multifractal  measures  resulting  from  multiplicative  processes  it  can 
be  shown  that  the  moments  of  the  normalized  measure  scales  according  to 


/  r  \  (q-l)Dq+d 

\L) 


(6.5) 


where  Dq  characterizes  the  scahng  of  the  qth  moment  (Meneveau  and  Sreenivasan,  1991). 
If  the  expected  value  in  6.5  is  obtained  from  the  measure  in  non-overlapping  segments  of 
size  r,  then  6.5  can  be  written 


/  r  \  {q-1)Dq 

\LJ 


(6.6) 


where  the  sum  is  taken  over  all  segments  of  length  r. 

Equation  6.6  can  be  used  to  estimate  Dq  by  raising  both  sides  to  the  1  /(q  -  1)  power; 


plotting 


yields  a  straight  line  with  a  slope  Dq. 

Regions  of  high  density  contribute  preferentially  to  moments  with  positive  q ,  and  re¬ 
gions  of  low  density  to  moments  of  negative  q.  As  |q|  increases,  moments  are  increasingly 
determined  by  the  extreme  behavior  of  the  measure,  by  the  highest  and  lowest  intensities, 
for  q  positive  and  negative  respectively.  The  moments  scale  with  r  as  determined  by  the 
exponents  Dq.  These  exponents  are  independent  of  the  scale  r  but  differ  for  moments  of 
different  order  q.  This  variation  is  characteristic  of  multifractal  measures;  exact  fractals,  by 
contrast,  have  identical  exponents  Dq  for  all  moments. 

The  two  families  of  exponents  we  have  presented  above  are  related:  from  the  curve  Dq 
one  can  obtain  the  singularity  spectrum  f(a),  and  vice-versa.  Each  order  q  provides  a  single 
(a,  /(a))  pair.  Define  r(q)  =  (q  -  1  )Dq;  then 


“<«>  =  ^T(?) 

/(<*(?))  =  qa(q)-T(q) 


For  a  derivation  of  6.7  see  Frisch  and  Parisi,  1985  or  Meneveau  and  Sreenivasan,  1991). 
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This  indirect  method  of  obtaining  the  singularity  spectrum  is  known  as  the  method  of 
moments.  We  will  use  it  below  to  explore  the  multifractal  structure  of  zooplankton  biomass 
variability.  (For  a  review  and  discussion  of  other  methods  see  Evertz  and  Mandelbrot, 
1992). 

6.3.2  Variables  analyzed 

The  biomass  time  series  6(f)  itself  is  at  best  only  weakly  multifractal  (more  details  below). 
Thus  we  also  analyzed  two  different  measures  of  the  variability  of  biomass: 


Si(t)  =  (b(t)  -  b{t  -  l))2 
S2(t)  =  (6(f) -<6))2 


where  brackets  denote  the  mean  value. 

These  quantities  estimate  total  variability  in  different  ways.  S2  is  the  familiar  squared 
deviation  from  the  mean;  it’s  expectation  is  the  variance  of  6,  and  measures  variability 
without  regard  to  temporal  autocorrelation.  The  expectation  of  Si  is  the  mean  square 
successive  difference,  which  measures  local  variability  in  consecutive  biomass  values.  It  is 
sensitive  to  autocorrelation  in  the  sequence:  if  the  series  is  positively  autocorrelated,  S j  will 
be  small,  and  vice-versa.  In  an  uncorrelated  random  series,  the  expectation  of  Si  is  twice 
the  variance  (von  Neumann  et  al.,  1941).  The  ratio  of  Si  to  S2  can  be  used  as  a  statistical 
test  for  a  first  order  Markov  process  against  the  alternative  of  random  variation. 

Figures  6-4  and  6-5  show  these  quantities  for  time  series  C  (Figure  6-1,  c).  To  in¬ 
vestigate  how  this  total  variability  is  organized  in  time,  we  subdivide  the  time  axis  into 
non-overlapping  intervals  of  length  r  and  compute  for  each  interval  the  following  normal¬ 
ized  measures 


Vi(t) 

V2(t) 


Sr  Sl(t) 
El  Si(t) 
E r_M) 
El  Suit) 


(6.8) 


where  J2r  denotes  the  sum  over  all  t  belonging  to  the  interval  of  length  r  centered  at  t, 
and  YIl  denotes  the  sum  over  the  whole  time  series.  Thus,  Vi  and  V2  are  normalized  sums 
of  squares  giving  the  proportion  of  the  total  variability  contributed  by  different  intervals  of 
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time. 


The  corresponding  normalized  densities  are 

-(t) = ffiH = (6-9) 

and 

"(t)  ■  = >>  (wo) 
The  following  argument  shows  that  tq  and  u2  are  in  fact  normalized  local  variances.  If 
N  measurements  occur  in  the  time  interval  L ,  and  n  in  each  interval  of  length  r,  then  L/r 
equals  N/n.  (Although  this  equality  appears  trivial  in  the  case  of  hourly  measurements,  it 
holds  for  any  frequency  of  sampling).  Thus  6.9  and  6.10  can  be  written  as 


»i(t)  = 


IZrSljt) 

F  El  Si  CO 


and 


V2(t)  =  -f 


k  Er  S2(0 


fElS2(0 


The  numerator  of  each  of  these  terms  is  an  average  squared  deviation;  i.e.,  a  variance,  within 
an  interval.  The  denominator  is  the  same  quantity  calculated  for  the  whole  data  set. 
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Figure  6-4:  Squared  first  differences  obtained  from  the  biomass  data  6(f)  in  time  series  C. 
(A)  4096  hours.  (B)  The  first  2048  hours  shown  separately  for  better  detail. 
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Figure  6-5:  Squared  differences  from  the  mean  for  the  biomass  data  b(t)  in  time  series  C. 
(A)  4096  hours.  (B)  The  first  2048  hours  shown  separately  for  better  detail. 

6.4  Analysis  of  zooplankton  variability 

We  begin  by  analyzing  the  longer  time  series  (series  C  in  Figure  6-1),  showing  that  both  V\ 
and  V2  are  multifractal  over  a  large  range  of  scales.  We  repeat  the  analysis  on  the  shorter 
time  series  (A  and  B)  to  investigate  the  generality  of  this  result. 

6.4.1  Scaling  of  moments 

We  consider  first  the  scaling  of  the  moments  for  V\  and  V2  (see  Figures  6-4  and  6-5).  The 
time  axis  is  subdivided  into  disjoint  intervals  of  length  r;  =  2*  (i  =  1, . . .,  11),  for  a  total 
length  L  =  4096  hours  (out  of  the  4099  hours  of  the  original  time  series).  For  each  scale 
n,  there  are  A;  =  L/21  intervals  over  which  to  compute  the  normalized  sums  of  squares  Vj 
and  V2.  The  sums  of  squares  in  interval  j  are  V\{j)  and  F2(j). 

If  equation  6.6  holds,  then  a  log-log  plot  of  the  (q  -  l)th  root  of  the  qth  moment  of  V\ 
or  V2  vs.  the  interval  length,  i.e.,  of 

E(u(;))' 

_j=i 

will  yield  a  straight  line  with  a  slope  Dq,  for  each  q.  The  same  will  be  true  for  V2. 

To  simplify  the  notation,  let  Pi(q,  r)  —  and  ^2(9,  r)  =  Sj(V2(j))l?.  Then 
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we  estimate  Dq  from  the  slope  of  a  plot  of  (log P\(q,r))/(q-  1)  vs.  log (r/L)  and  similarly 
for  P2.  These  plots  are  shown  in  Figures  6-6  and  6-7  for  some  representative  values  of  q  in 
[—3,  +3]  . 


log  r/L 


log  r/L 


Figure  6-6:  Scaling  of  moments  for  measure  VI.  Log  Pi(q,  r)?-1  vs.  log  r/L  for  represen¬ 
tative  values  of  q:  (a):  q  =  —3,  (b):  q  =  —  1,  (c):  q  =  +1.25,  and  (d):  q  =  +3. 
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Figure  6-7:  Scaling  of  moments  for  measure  V2.  Log  P2(q,r)<i-1  vs.  log  r/L  for  represen¬ 
tative  values  of  q:  (a):  q  —  —3,  (b):  q  =  -1,  (c):  q  =  +1.25,  and  (d):  q  =  +3. 
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The  slope  Dq  must  be  estimated  only  over  the  range  of  scale  values  for  which  the 
curve  is  linear;  this  range  is  known  as  the  scaling  region.  In  Figures  6-6  and  6-7  (  and 
equivalent  plots  for  the  other  q  values),  the  smallest  scale  at  which  the  curves  display 
linearity  appears  to  lie  between  23  and  24,  i.e.  eight  to  sixteen  hours.  This  limit  may  result 
from  the  processes  generating  the  data,  from  noise  in  the  measurements  or  from  problems 
with  moment  convergence  at  high  and  low  q  values.  Noise  is  known  to  produce  curvature 
at  small  scales  for  the  most  negative  q  values  (Meneveau  and  Sreenivasan,  1991),  and  our 
data  are  more  linear  the  closer  q  is  to  zero  (see  Figures  6-8  and  6-9). 

We  chose  r  =  23  as  the  lower  limit  and  r  =  211  as  the  upper  limit  of  the  scaling  region, 
and  fit  straight  lines  to  the  data  by  least-squares.  The  lines  fit  the  data  well,  with  coefficients 
of  determination  R 2  =  0.997  ( q  =  —3)  and  R 2  =  0.971  (q  =  +3)  for  Vj,  and  R2  =  0.985 
(q  =  -3)  and  R2  =  0.99  (q  =  +3)  for  V2.  Figures  6-8  and  6-9  show  these  log-log  plots  in 
the  scaling  region  for  selected  values  of  q.  The  slopes  of  these  lines  are  the  exponents  Dq. 


Figure  6-8:  Scaling  region  for  moments  of  FI.  The  lines  are  the  least-square  fit  trough  the 
data  points.  From  bottom  to  top  q  values  are  -3,  -2,  -1,  0,  +1.25,  +2,  +3. 
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Figure  6-9:  Scaling  region  for  moments  of  V2.  The  lines  are  the  least-square  fit  trough  the 
data  points.  From  bottom  to  top  q  values  are  -3,  -2,  -1,  0,  +1.25,  +2,  +3. 

We  plot  in  Figure  6-10  (a,b)  the  slopes  Dq  as  a  function  of  q  for  both  V\  and  V2.  The 
variation  of  Dq  with  q  is  characteristic  of  multifractal  structures.  This  variation,  coupled 
to  the  linear  scaling  of  the  moments  in  a  large  range  of  temporal  scales,  provides  evidence 
for  the  multifractal  structure  of  hi  and  V2. 

6.4.2  The  singularity  spectrum 

From  the  scaling  of  the  moments,  the  pairs  (a,  /(a))  are  computed  via  the  transformations 
in  6.7.  The  derivative  of  (q  -  1  )Dq  was  estimated  by  centered  differences  for  q  values  at 
intervals  of  0.25  in  [-3, +3].  The  resulting  /(a)  curves  are  shown  in  Figure  6-11  (a,b)  for 
Vi  and  V2  respectively.  These  curves  display  a  parabolic  shape  characteristic  of  multifractal 
measures.  The  maxima  for  /(a)  corresponds  to  q  =  0  and  equals  the  Euclidean  dimension 
of  the  measure  support  (here,  equal  to  one). 
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Figure  6-12:  Comparison  of  the  curve  f(a)  for  V2  in  three  different  data  sets:  (+),  time 
series  A,  (o),  time  series  B,  and  (*),  time  series  C. 

Just  as  the  variation  in  Dq  values  reveals  the  inhomogeneity  of  Vf  and  Vi  along  the  time 
axis,  so  does  the  spread  in  a  values  around  the  maximum  of  /(a).  Both  are  characteristic 
of  multifractal  measures.  We  have  compared  these  results  to  those  for  other  intermittent 
quantities  described  as  multifractals,  and  to  data  obtained  numerically  from  the  binomial 
process  (Meneveau  and  Sreenivasan,  1991;  Prasad  et  al. ,  1988).  The  spread  observed  here 
lies  well  within  the  range  of  these  other  studies.  We  also  compared  the  results  with  two 
equally  long  data  sets  known  not  to  be  multifractal:  white  noise  and  a  time  series  of 
velocity  in  a  turbulent  flow  (provided  by  K.R.  Sreenivasan).  For  q  in  [—3,  +3],  the  maximum 
difference  in  Dq  values  was  0.02  for  white  noise,  and  0.09  for  the  velocity  data,  well  below 
the  values  of  0.73  and  0.83  for  V\  and  Vi- 

We  repeated  this  analysis  for  Vi  with  the  two  other  time  series  (A,B).  These  data  sets 
are  shorter  (2732  and  2735  data  points  respectively);  we  used  only  the  last  2048  (or  211) 
points.  The  scaling  of  the  moments  holds,  and  the  exponents  Dq  were  obtained  for  the  same 
scaling  region  as  above.  Figure  6-12  compares  the  resulting  f(a)  curves  with  the  curve  for 
time  series  C.  The  three  curves  are  very  similar,  particularly  for  a  <  1  (q  >  0). 

Finally,  the  spikiness  of  the  biomass  data  suggests  the  possibility  of  the  biomass  itself 
being  multifractal.  We  have  studied  this  question  and  the  results,  not  presented  here, 
indicate  that  biomass  does  not  present  a  convincing  multifractal  structure.  The  scaling  of 
the  moments  presents  a  large  scaling  region  but  the  variation  in  the  resulting  exponents 
(0.16)  is  small  compared  to  other  data  sets  described  as  multifractals. 
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6.4.3  Robustness  of  conclusions 

Could  our  findings  be  an  artifact  of  infrequent  large  peaks  located  randomly  in  a  data 
set  of  length  LI  Multifractal  processes  include  a  much  richer  structure  than  this  (e.g., 
the  binomial  process  in  Section  6.3.1).  As  a  test  for  such  artifacts,  we  generated  data  by 
randomly  permuting  the  time  series  (b(t)  -  b(t  -  l))2  and  (b(t)  -  ( b ))2.  Any  multifractal 
structure  in  the  real  data  should  be  destroyed  by  this  procedure.  In  fact,  the  permuted 
time  series  displayed  a  much  smaller  scaling  range  (only  4  orders  of  magnitude,  from  27  to 
210)  and  a  reduced  spread  of  Dq  values  over  this  range.  For  30  permuted  time  series  of  Vi, 
and  q  in  [-3,  +3],  the  mean  range  in  Dq  was  0.275,  and  the  maximum  range  was  0.39.  The 
corresponding  values  for  V2  were  0.18  and  0.25. 

Could  our  results  be  due  to  the  limited  length  of  the  time  series?  L  =  4096  seems  large 
for  an  ecological  time  series,  but  is  short  compared  to  data  sets  that  have  been  described  as 
multifractals  in  other  fields.  The  length  of  the  data  set  limits  the  range  of  orders  q  that  can 
be  considered.  Moments  with  high  positive  or  low  negative  q  are  determined  by  the  extreme 
behavior  of  the  data.  Since  such  extreme  behavior  appears  infrequently  in  an  intermittent 
signal,  convergence  of  statistical  estimates  of  moments  with  large  \q\  requires  a  long  record. 
This  convergence  is  necessary  to  obtain  reliable  estimates  of  the  exponents  Dq. 

Equation  6.5  shows  that  the  relevant  quantities  in  calculating  the  exponents  Dq  are  the 
logarithms  of  the  moments  divided  by  (q  -  1).  To  explore  convergence  of  these  quantities, 
we  studied  the  behavior  of  the  moments  of  the  densities  xq  and  v2  as  a  function  of  time 
series  length,  following  the  approach  of  Meneveau  and  Sreenivasan  (1991).  We  consider 
moments  of  vi  rather  than  of  Vi  {i  =  1,2)  because  as  L  becomes  large  the  former  tend  to 
a  constant  value  while  the  latter  decrease  with  L  (Meneveau  and  Sreenivasan,  1991).  The 
two  are  related  by  ((v,-)9)  =  (V?)(L/r)q. 

The  behavior  of  the  moments  with  length  L  for  the  two  extreme  values  of  q  (-3  and 
+3)  is  shown  in  Figures  6-13  and  6-14.  The  different  curves  correspond  to  different  scales  r. 
Although  the  moments  have  not  converged  for  L  =  4096,  the  differences  between  the  curves 
become  fairly  constant  as  L  increases.  Hence,  the  ratio  between  moments  at  two  different 
scales  tends  to  a  constant  independent  of  time  series  length.  It  is  exactly  these  ratios  that 
are  used  to  obtain  the  slopes  Dq.  Their  convergence  thus  supports  our  conclusion  that  the 
time  series  is  long  enough  for  the  interval  of  q  values  used  here. 
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Figure  6-13:  The  logarithm  of  the  moments  (vl9)  divided  by  (q  -  1)  as  a  function  of  time 
series  length  L  for  two  different  values  of  q:  (A)  q  =  -3,  (B)  q  —  +3.  The  different  curves 
in  each  plot  correspond  to  different  scales  r  =  2!.  (In  (A),  i  varies  from  3  to  10,  in  (B), 
from  10  to  3  for  curves  from  bottom  to  top). 
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Figure  6-14:  The  logarithm  of  the  moments  ( v2q )  divided  by  (q  —  1)  as  a  function  of  time 
series  length  L  for  two  different  values  of  q:  (A)  q  =  -3,  (B)  q  =  +3.  The  different  curves 
in  each  plot  correspond  to  different  scales  r  =  2‘.  (In  (A),  i  varies  from  3  to  10,  in  (B), 
from  10  to  3  for  curves  from  bottom  to  top). 


6.5  Discussion 


Our  analyses  show,  we  think  convincingly,  that  zooplankton  biomass  variability  is  highly 
intermittent,  and  can  be  characterized  as  a  multifractal.  Multifractal  analysis  is  a  recent 
development,  and  this  is  its  first  application  to  zooplankton.  Because  of  this,  assessing  the 
significance  of  our  findings  raises  more  questions  than  it  answers.  We  consider  some  of  those 
here. 

Like  most  oceanographic  data,  the  data  analyzed  here  contain  both  spatial  and  temporal 
components,  because  biomass  was  measured  at  a  fixed  location  in  a  moving  water  mass. 
More  work  is  needed  to  determine  the  relative  contribution  of  spatial  heterogeneity  and 
temporal  dynamics  to  the  intermittency.  In  a  different  analysis  of  the  same  data  set,  Flagg  et 
al.  (1994)  could  not  explain  zooplankton  fluctuations  by  advection  and  spatial  heterogeneity 
alone,  which  suggests  that  both  might  be  important.  They  report  that  the  acoustic  data 
were  collected  near  a  front  separating  offshore  and  inshore  water  masses,  and  that  the 
position  of  this  front  was  affected  by  large  scale  cross-shelf  movements  (approximately  40 
Km  and  a  month  duration),  wind-driven  advective  events  with  a  typical  time  scale  of  days, 
and  tidal  advection  (Flagg  et  al .,  1994;  Wirick,  1994).  There  was  limited  coherence  between 
this  physical  variability  and  zooplankton  biomass,  perhaps  because  of  the  interaction  of 
the  physics  with  biological  processes,  including  predator-prey  interactions,  migration  and 
aggregation  (Ascioti  et  al .,  1993;  Flagg  et  al,  1994).  These  physical  and  biological  processes 
operate  over  a  range  of  scales  comparable  to  the  scaling  range  we  have  found. 

We  conjecture  that  intermittent,  multifractal  patterns  will  eventually  be  found  in  both 
spatial  variability  and  the  temporal  variability  of  spatial  averages  in  zooplankton  data.  It 
is  only  fair  to  ask,  “so  what?”;  how  might  multifractal  data  analysis  be  applied  to  problems 
in  plankton  ecology.  Four  applications  spring  to  mind:  designing  sampling  schemes,  char¬ 
acterizing  the  environment  encountered  by  planktonic  organisms,  comparative  studies,  and 
inferring  processes. 

•  Intermittency  has  important  implications  for  sampling  (Bohle-Carbonell,  1992).  Un¬ 
dersampling  an  intermittent  signal  is  particularly  problematic;  it  leads  to  biased  esti¬ 
mates  of  means  and  confidence  limits  (Baker  and  Gibson,  1987).  Simple  multiplicative 
cascade  models  could  provide  a  useful  tool  to  simulate  intermittent  data  with  a  given 
multifractal  structure,  for  evaluation  of  potential  sampling  schemes. 


159 


•  The  environment  is  sampled  not  only  by  oceanographers  but,  in  a  sense,  by  the  or¬ 
ganisms  that  live  in  the  ocean.  The  small-scale  variability  experienced  by  individual 
planktonic  organisms  may  have  important  implications  for  foraging,  behavior,  growth, 
and  population  dynamics  (Davis  and  Flierl,  1991;  Goldman,  1988;  Rothschild,  1992). 
If  the  environment,  either  physical  or  biological,  is  multifractal,  the  relations  between 
variability  and  scale  will  be  richer  and  more  structured  than  otherwise.  The  conse¬ 
quences  of  this  richness  for  organisms  living  in  this  environment  remain  to  be  explored. 

•  The  initial  uses  of  multifractal  analysis,  like  any  other  descriptive  statistic,  will  be 
exploratory  and  comparative.  It  will  require  more  examples  to  discover  how  to  use 
the  analysis  to  detect  interesting  differences  between  taxa,  habitats,  locations,  sea¬ 
sons,  etc.  Spectral  analysis  has  long  been  used  in  this  way.  Multifractals  provide 
a  complementary  approach.  While  spectral  analysis  investigates  the  relative  contri¬ 
bution  of  different  scales  to  total  variance,  multifractals  reveal  the  organization  or 
structure  of  variability  in  space  or  time.  The  relationship  of  simple  fractals  to  the 
variance  power  spectrum  is  known  (Rothrock  and  Thorndike,  1980),  that  of  multi¬ 
fractals  is  not.  Eventually,  one  would  hope  that  a  relation  between  multifractals  and 
higher-order  spectral  analysis  would  be  established. 

Comparisons  can  be  based  directly  on  the  curves  /(a)  or  D{q),  or  on  other  summary 
indices  calculated  from  such  curves  (Prasad  et  al.,  1988).  One  such  index  is  amin,  the 
minimum  value  of  a,  for  which  /(a)  =  0.  This  value  of  the  local  exponent  a  measures 
the  highest  degree  of  singularity  in  the  data  (i.e.  the  most  spiky  behavior;  see  Eq 
6.1).  By  comparing  am,n  one  can  determine  which  data  set  contains  the  strongest 
singularities  (see  Prasad  et  al.  1988)  for  an  example).  In  practice,  am{n  is  calculated 
by  extrapolating  the  curve  /(a),  since  the  length  of  the  data  set  limits  the  range  of  a 
values  one  can  obtain. 

Another  potentially  useful  index  is  the  intermittency  exponent  fi,  which  quantifies 
the  degree  of  intermittency  of  the  data.  The  intermittency  exponent  is  determined 
by  the  spread  of  a  around  its  mean.  Recall  that  r(g)  =  (q  —  1  )Dq  and  let  fi  = 
— d?T(q) / dq2\q=o-  Then  the  variance  of  a  is  given  by 

al  =  v/HL/r) 
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(Meneveau  and  Sreenivasan,  1991).  The  larger  the  value  of  /x,  the  larger  the  variance 
of  a  and  the  higher  the  degree  of  intermittency  of  the  data.  This  can  also  be  seen  by 
the  relation  of  /z  to  the  variance  of  the  measure  In (Mt/Ml) 

°Mr  =  VHL/r) 

The  intermittency  exponent  [i  was  introduced  by  Kolmogorov  (1962)  in  studies  of  the 
energy  dissipation  rate  in  turbulence.  For  the  longer  time  series  analyzed  here,  the 
intermittency  exponent  of  Vj  (fi  =  0.3)  is  higher  than  that  of  V<i  ( fj,  =  0.17)  (both 
computed  by  centered  differences  from  the  curve  r(q)).  The  former  is  comparable  to 
values  in  the  literature  for  the  intermittency  exponents  of  energy  and  scalar  dissipation 
rates  in  turbulent  flows  (Prasad  et  al.  ,  1988).  In  the  ocean,  strong  intermittency  has 
been  found  in  physical  quantities  such  as  the  dissipation  rates  of  turbulent  velocity  and 
temperature  fluctuations  (Baker  and  Gibson,  1987).  The  implications  of  such  physical 
intermittency  for  the  distribution  of  biological  variables  remain  to  be  explored. 

The  multifractal  formalism  has  been  extended  to  more  than  one  variable  to  describe 
the  degree  of  correlation  in  intermittent  quantities  in  turbulence  (Meneveau  et  al ., 
1989).  This  suggest  that  multifractal  analyses  could  also  be  used  to  compare  the 
distributions  of  plankton  variability  and  those  of  passive  scalars  and  environmental 
variables. 

•  One  would  hope  that,  eventually,  multifractal  statistics  would  help  identify  the  pro¬ 
cesses  producing  the  pattern.  At  the  present,  this  is  not  possible.  It  is  known  that 
multifractals  can  be  produced  by  multiplicative  processes,  and  that  they  appear  in 
the  measures  of  trajectories  on  strange  attractors.  This  does  not,  however,  mean 
that  a  multifractal  zooplankton  pattern  is  produced  by  a  multiplicative  process  or 
as  a  strange  attractor.  We  simply  do  not  know  enough  about  the  possibilities  for 
producing  multifractals.  We  expect  that  the  connection  between  pattern  and  pro¬ 
cess  for  multifractal  variability  in  the  plankton  will  develop  along  a  similar  path  to 
that  of  spectral  analysis.  The  initial  uses  of  spectral  analysis  were  purely  descriptive 
(Platt,  1972;  Platt  et  al.,  1970).  That  use  was  followed  by  a  connection  of  spectral 
analysis  to  phenomenological  models  (cascade  models  of  how  variance  transfers  from 
larger  to  smaller  scales;  Denman  and  Platt,  1976),  and  only  later  by  a  connection  to 
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mechanistic  models. 


However,  it  may  be  possible  to  identify  spatial  or  temporal  scales  over  which  the 
processes  must  be  different,  by  identifying  limits  to  the  scaling  region.  This  occurs 
in  studies  of  turbulence,  where  the  scaling  regime  is  different  above  and  below  the 
Kolmogorov  scale  for  the  dissipation  rate  of  a  passive  scalar.  This  change  reflects  the 
different  dominant  physical  processes  operating  in  these  two  regions  (Sreenivasan  and 
Prasad,  1989).  A  well  known  biological  example  is  the  change  in  the  spectral  exponent 
that  occurs  at  the  so-called  ‘Platt  knee’  in  the  power  spectrum  of  phytoplankton 
spatial  data.  This  change  has  been  related  to  a  switch  from  the  influence  of  the 
physical  factor  of  turbulence  to  the  biological  factor  of  reproduction  (Denman  and 
Platt,  1976). 

These  applications  of  multifractals  await  further  investigations  of  other  zooplankton 

data  sets.  This  paper  has  presented  evidence  for  the  potential  of  multifractals  to  become 

an  important  descriptive  tool  in  zooplankton  ecology. 
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Afterthoughts 


Comunque  il  signor  Palomar  non  si  perde  d’animo  e  a  ogni  momento  crede 
d’esser  riuscito  a  vedere  tutto  quel  che  poteva  vedere  dal  suo  punto 
d’osservazione,  ma  poi  salta  fuori  sempre  qualcosa  di  cui  non  aveva  tenuto 
conto.  Se  non  fosse  per  questa  sua  impazienzia  di  raggiungere  un  risultato 
completo  e  definitivo  della  sua  operazione  visiva,  il  guardare  le  onde  sarebbe 

per  lui  un  esercizio  molto  riposante .  E  forse  potrebbe  essere  la  chiave  per 

padroneggiare  la  complessita  del  mondo  riducendolo  al  meccanismo  piu 
sernplice. 

....Il  signor  Palomar  s’allontana  lungo  la  spiaggia,  coi  nervi  tesi  com’era 
arrivato  e  ancor  piu  insicuro  di  tutto ? 

— Italo  Calvino.  Palomar. 

In  each  of  the  preceding  chapters,  I  have  presented  specific  conclusions  and  directions 
for  future  research.  Here,  I  return  to  the  main  theme  presented  in  the  introduction:  the 
mismatch  of  environmental  and  biological  scales  in  nonlinear  ecological  systems.  I  state 
some  main  conclusions  of  the  work  in  light  of  that  theme,  and  briefly  speculate  on  possible 
future  developments. 

This  work  has  demonstrated  two  novel  ways  by  which  a  scale  mismatch  would  occur: 

•  The  first  one  relates  to  population  structure  in  consumer- resource  interactions.  The 
distribution  of  the  consumer  population  in  life-history  stages  is  potentially  impor- 

2 My  own  translation:  ‘Anyhow  Mister  Palomar  does  not  lose  heart  and  at  every  moment  believes  he  has 
succeeded  at  seeing  everything  that  he  could  see  from  his  observation  point,  but  then  something  he  had  not 
considered  always  comes  up.  If  it  were  not  for  his  impatience  to  reach  a  complete  and  definite  result  of  his 
visual  operation,  looking  at  the  waves  could  be  a  very  relaxing  exercise... And  it  could  be  the  key  to  master 
the  complexity  of  the  world  by  reducing  it  to  the  simplest  mechanism.. .Mister  Palomar  goes  away  along  the 
shore,  with  his  nerves  tense  as  he  had  arrived  and  even  more  insecure  of  everything.’ 
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tant  to  consumer-resource  dynamics.  Specifically,  this  distribution  can  lead  to  cy¬ 
cles  of  the  consumer-resource  interaction  under  a  constant  resource  supply.  When 
these  endogenous  cycles  interact  with  the  frequencies  of  a  variable  resource  supply, 
the  temporal  scales  in  consumer-resource  dynamics  can  differ  from  those  in  the  forc¬ 
ing.  Unstructured  models  may  therefore  miss  an  essential  element  of  the  response  of 
consumer-resource  interactions  to  environmental  forcings. 

•  The  second  one  relates  to  the  spatial  coupling  of  predator-prey  cycles  in  heterogeneous 
space.  Such  a  coupling  can  lead  to  complex  spatio-temporal  dynamics  (i.e.  chaos  and 
quasiperiodicity)  of  the  predator  and  prey.  In  these  dynamic  regimes,  the  spatial 
patterns  of  the  populations  differ  in  scale  from  the  underlying  gradient.  Essential 
ingredients  for  such  a  scale  mismatch  are  the  local  endogenous  cycles  of  predator 
and  prey.  Complex  dynamics  results  from  the  weak  coupling  by  diffusion  of  cycles 
that  differ  in  space  because  of  the  underlying  gradient.  These  results  emphasize  that 
studies  of  the  physical  forcing  of  food  webs  may  be  extremely  sensitive  to  the  local 
biological  dynamics.  In  heterogeneous  environments,  the  coupling  of  local  limit  cycles 
may  lead  to  drastically  different  patterns  than  the  coupling  of  local  equilibria. 

Simple  models,  such  as  the  ones  studied  here,  identify  potential  scenarios  for  the  mismatch 
of  scales  in  nature.  I  anticipate  a  wealth  of  surprising  new  results  in  this  area,  with  both 
spatial  and/ or  stochastic  models  of  nonlinear  ecological  interactions.  The  epidemiological 
models  described  briefly  in  the  introduction,  already  give  us  a  glimpse  of  the  unexpected 
consequences  of  stochasticity  in  nonlinear  systems.  Much  more  will  probably  come  from 
studies  of  the  interplay  of  stochasticity  with  multiple  attractors  and/or  repellors.  Stochastic 
forcings  are  ubiquitous  in  nature,  they  represent  not  only  the  high  dimensional  fluctuations 
in  the  environment  but  also  the  effect  of  low  population  numbers  on  dynamics.  The  coex¬ 
istence  of  attractors  and  repellors  may  also  be  common  in  ecological  systems. 

Theoretical  results  on  forced  ecological  models  raise  an  important  empirical  question: 
how  can  we  identify  environmental  forcings  related  to  specific  biological  patterns  when 
scales  do  not  match?  Some  recent  methods  in  nonlinear  data  analysis  could  provide  a 
tool  to  approach  this  question.  These  methods  were  developed  by  Ellner  and  Turchin 
(1995)  to  compute  Lyapunov  exponents  from  ecological  time  series  that  are  both  noisy  and 
short.  They  are  based  on  the  idea  of  ‘reconstructing’  the  dynamics  of  a  multidimensional 
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system  from  time-delay  coordinates  of  a  single  variable  (see  Kot  et  al.,  1988,  for  ecological 
discussion).  Recently,  ‘reconstruction’  was  extended  to  stochastic  nonautonomous  systems 
(Casdagli,  1992).  Based  on  these  ideas,  the  methods  of  Ellner  and  Turchin  (1995)  fit,  to  a 
univariate  data  set  Nt,  variants  of  the  following  time-series  model 


Nt+L  = 

where  d  denotes  the  embedding  dimension  (how  far  in  the  past  one  must  look  for  an  ex¬ 
planation  of  current  changes  in  A),  and  X,  an  arbitrary  time  lag.  The  two  terms  et  and 
Et  are  a  stochastic  term,  representing  dynamic  noise,  and  a  periodic  function,  representing 
seasonality.  These  two  terms  make  the  fitted  model  nonautonomous,  they  estimate  the 
exogenous  components  of  the  system.  Without  these  terms,  the  function  /  represents  the 
endogenous  structure  of  the  system.  These  methods  could  provide  a  basis  to  estimate  the 
qualitative  behavior  of  endogenous  structures,  and  to  test  hypotheses  on  periodic  forcings 
at  specific  temporal  scales.  Identifying  the  type  of  endogenous  dynamics  is  a  first  step  to 
establish  the  potential  consequences  of  temporal  forcings. 

Casdagli,  M.  1992.  A  dynamical  systems  approach  to  modeling  input-output  systems.  In 
M.  Casdagli  and  S.  Eubank,  eds.,  Nonlinear  modeling  and  forecasting.  SFI  Studies  in 
the  Sciences  of  Complexity  Proc.  Vol.  XII.  Addison- Wesley,  New  York. 

Ellner,  S.  and  P.  Turchin.  1995.  Chaos  in  a  ’noisy’  world:  new  methods  and  evidence  from 
time  series  analysis,  in  press. 

Kot,  M.,  W.M.  Schaffer,  G.L.  Truty,  D.J.  Graser,  and  L.F.  Olsen.  1988.  Changing  criteria 
for  imposing  order.  Ecol.  Model.  43:  75-110. 

Sirovich,  L.  1987.  Turbulence  and  the  dynamics  of  coherent  structures.  Part  I:  coherent 
structures.  Quart.  Appl.  Math.  XLV(3):  561-571. 
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Appendix 


Relationships  among  measures  of  characteristic  scale. 

Although  different  measures  of  characteristic  scale  can  be  found  in  the  ecological  lit¬ 
erature,  the  relations  among  these  quantities  are  missing.  I  present  below  four  common 
measures  and  show  how  they  all  relate  through  the  autocorrelation  function.  These  re¬ 
lationships  illustrate  that  all  four  quantities  are  essentially  measuring  the  same  intuitive 
concept  of  scale:  the  distance  (or  time)  one  has  to  travel  (or  wait)  to  see  a  significant 
change  in  the  variable  of  interest  (Powell,  1989). 

The  four  measures  of  characteristic  scale  are: 


I  Correlation  length. 

The  value  of  the  lag  at  which  the  autocorrelation  function  first  crosses  zero.  The 
theoretical  autocorrelation  function  of  a  stationary  stochastic  process  Y(x)  is  given 


by 


m  _  E{lY(x)~  fAiY(x  +  0-/*]} 

P[  >  E{[Y(x )  -  p)*} 


(6.11) 


where  E  denotes  expected  value,  /,  a  lag,  and  p  =  E{Y).  In  words,  p(l )  is  the  ratio 
between  the  covariance  of  values  separated  by  a  lag  l  and  the  variance  of  Y .  Clearly, 
p(0)  =  1,  p{l)  initially  decreases  with  l,  and  p[l)  varies  between  1  and  —1.  (For 
empirical  estimates  of  the  autocorrelation  function  see  Chatfield  (1989)  or  Chapter 
5). 


The  correlation  length  is  given  by  the  value  l  at  which  p(l)  =  0.  Note  that  this 
definition  measures  a  ‘significant  change’  in  the  quantity  of  interest,  by  a  significant 
decrease  in  the  autocorrelation  function. 
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II  Period  or  wavelength  associated  with  a  peak  in  the  power  spectrum. 

This  definition  is  widely  used  in  biological  oceanography.  Because  the  power  spectrum 
gives  the  distribution  of  variance  as  a  function  of  frequency  or  wavenumber,  a  peak 
in  the  spectrum  indicates  a  large  proportion  of  the  variance  at  the  associated  period 
or  wavelength. 


Ill  The  microscale. 

Powell  (1989)  introduced  into  ecology  a  measure  of  scale  known  in  turbulence  as  the 
microscale.  Is  is  given  by 


Lr 


E{(Y(x)-rf} 


(6.12) 


The  scale  Lm  is  obtained  from  the  ratio  between  the  variance  of  Y  and  its  mean 
squared  derivative.  For  empirical  data,  the  mean  squared  derivative  can  be  evaluated 
as  the  mean  of  the  squared  first  differences. 


IV  The  integral  scale. 

Another  proposed  measure  of  scale  is  the  integral  of  the  autocorrelation  function, 

Li  =  jp(l)dX.  (6.13) 

For  an  application  to  biological  oceanography,  see  Mackas  (1984).  At  first,  both  the 
microscale  and  the  integral  scale  are  difficult  to  interpret.  Their  relations  to  the 
autocorrelation  function  will  clarify  how  they  provide  a  measure  of  scale. 

The  relationship  between  measures  (I)  and  (II)  is  best  known.  In  fact,  the  power 
spectrum  is  the  cosine  Fourier  transform  of  the  autocorrelation  function  (see  Platt  and 
Denman  1975  for  a  presentation  in  the  ecological  literature  of  both  functions  and  their 
relationship).  Periodicity  in  the  data  leads  to  wavelike  peaks  in  the  correlogram,  and  the 
correlogram  crosses  zero  at  a  lag  equal  to  a  fourth  of  the  period. 

The  relationship  between  measures  (I)  and  (III)  is  as  follows.  If  the  autocorrelation 
function  is  approximated  at  the  origin  by  a  parabola,  the  scale  Lm  gives  the  value  of  /  for 
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which  the  parabola  decreases  from  1  to  1/2.  Hence,  a  ‘significant  change’  in  the  quantity 
of  interest  is  measured  by  a  significant  decrease  in  the  autocorrelation  function,  or  more 
exactly,  in  the  parabola  that  approximates  the  autocorrelation  function  at  small  lags.  This 
is  shown  below  for  a  variable  Y(x )  whose  mean  ji  =  0.  This  choice  eliminates  cumbersome 
notation;  the  argument,  however,  holds  for  any  stationary  variable  with  a  mean  other  than 
zero.  With  fi  =  0,  the  correlation  function  becomes 


E{Y(x  +  l)Y(x)} 
P[  j  E{Y\x)} 

It  can  be  written  as 


(6.14) 


E{Y\x  +  /)}  +  E{Y2(x)}  -  E{[Y(x  +  /)  -  Y{x)}2} 
2  E{Y2(x)} 


(6.15) 


But,  because  Y(x)  is  stationary,  equation  6.15  simplifies  to 

,  =  2 E{Y\X)}  -  -£{[r(*  +  l)  -  Y(x)f} 

2E{Y\x )} 

’“2  WWf}  ' 


If  l  becomes  arbitrarily  small,  then 


limp(/)  = 


l2 

2  E{Y\x )} 


2/4/ 


Thus,  for  small  /,  the  correlation  function  can  be  approximated  by  the  parabola  y(l)  = 
1  -  l2!2L2m.  This  parabola  crosses  zero  at  /  =  y/2Lm  and  decreases  from  1  to  1/2  for 
l  =  Lm.  This  completes  the  argument  relating  Lm  to  the  autocorrelation  function. 

Finally,  the  integral  scale  and  the  autocorrelation  function  are  trivially  related  by  defi¬ 
nition.  The  interpretation  of  the  integral  scale  can  be  clarified  with  the  following  argument 
borrowed  from  McComb  (1990).  Approximate  the  autocorrelation  function  by  a  decreasing 
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exponential  function  p(l )  ~  exp(-6/).  Then, 

Lj  —  J  exp(—bl)dl  =  -.  (6.16) 

But,  p(l/b)  -  1/e.  Thus,  the  integral  scale  corresponds  to  the  lag  for  which  the  auto¬ 
correlation  function  equals  1/e,  provided  the  autocorrelation  function  is  well  approximated 
by  a  decreasing  exponential.  Therefore,  the  integral  scale  measures  a  ‘significant  change’  in 
the  variable  of  interest  by  a  decrease  in  the  autocorrelation  function  from  1  to  1/e. 

In  summary,  all  four  measures  of  scale  relate  to  the  autocorrelation  function.  Roughly 
speaking,  they  tell  us  how  far  to  travel  (or  how  long  to  wait)  for  the  data  to  become 
uncorrelated  with  itself. 

Chatfield,  C.  1975.  The  analysis  of  time  series,  an  introduction.  Chapman  and  Hall,  New 
York. 

Mackas,  D.L.  1984.  Spatial  autocorrelation  of  plankton  community  composition  in  a  conti¬ 
nental  shelf  ecosystem.  Limnology  and  oceanography  29(3):  451-471. 

McComb,  W.D.  1990.  The  Physics  of  Fluid  Turbulence.  Clarendon  Press.  Oxford. 

Platt,  T.  and  Denman,  K.L.  (1975)  Spectral  analysis  in  ecology.  Ann.  Rev.  Ecol.  Syst.  6: 
189-210. 

Powell,  T.  M.  1989.  Physical  and  biological  scales  of  variability  in  lakes,  estuaries,  and 
the  coastal  ocean.  In  J.  Roughgarden,  R.M.  May  and  S.  Levin,  eds.,  Perspectives  in 
ecological  theory.  Princeton  Univ.  Press,  Princeton,  New  Jersey. 
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